{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating Thermodynamics Observables with a quantum computer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "\n",
    "from qiskit.utils import QuantumInstance\n",
    "from qiskit import Aer\n",
    "from qiskit.algorithms import NumPyMinimumEigensolver, VQE\n",
    "\n",
    "from qiskit_nature.drivers import UnitsType, Molecule\n",
    "from qiskit_nature.drivers.second_quantization import (\n",
    "    ElectronicStructureDriverType,\n",
    "    ElectronicStructureMoleculeDriver,\n",
    ")\n",
    "from qiskit_nature.problems.second_quantization import ElectronicStructureProblem\n",
    "from qiskit_nature.converters.second_quantization import QubitConverter\n",
    "from qiskit_nature.mappers.second_quantization import JordanWignerMapper\n",
    "from qiskit_nature.algorithms import GroundStateEigensolver\n",
    "import qiskit_nature.constants as const\n",
    "from qiskit_nature.algorithms.pes_samplers import BOPESSampler, EnergySurface1DSpline\n",
    "from thermodynamics_utils.thermodynamics import constant_volume_heat_capacity\n",
    "from thermodynamics_utils.vibrational_structure_fd import VibrationalStructure1DFD\n",
    "from thermodynamics_utils.partition_function import DiatomicPartitionFunction\n",
    "from thermodynamics_utils.thermodynamics import Thermodynamics\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.simplefilter(\"ignore\", np.RankWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A preliminary draft with more information related to this tutorial can be found in preprint: Stober et al, arXiv 2003.02303 (2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation of the Born Oppenheimer Potential Energy Surface (BOPES) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute thermodynamic observables we begin with single point energy calculation which calculates the wavefunction and charge density and therefore  the energy of a particular arrangement of nuclei. Here we  compute the Born-Oppenheimer potential energy surface of a hydrogen molecule, as an example, which is simply the electronic energy as a function of bond length. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubit_converter = QubitConverter(mapper=JordanWignerMapper())\n",
    "quantum_instance = QuantumInstance(backend=Aer.get_backend(\"aer_simulator_statevector\"))\n",
    "solver = VQE(quantum_instance=quantum_instance)\n",
    "\n",
    "me_gss = GroundStateEigensolver(qubit_converter, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))\n",
    "mol = Molecule(\n",
    "    geometry=[(\"H\", [0.0, 0.0, 0.0]), (\"H\", [0.0, 0.0, 0.2])],\n",
    "    degrees_of_freedom=[stretch1],\n",
    "    masses=[1.6735328e-27, 1.6735328e-27],\n",
    ")\n",
    "\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = ElectronicStructureMoleculeDriver(mol, driver_type=ElectronicStructureDriverType.PYSCF)\n",
    "\n",
    "es_problem = ElectronicStructureProblem(driver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BOPES sampler testing\n",
    "bs = BOPESSampler(gss=me_gss, bootstrap=True)\n",
    "points = np.linspace(0.45, 5, 50)\n",
    "res = bs.sample(es_problem, points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = []\n",
    "bs_res_full = res.raw_results\n",
    "for point in points:\n",
    "    energy = bs_res_full[point].computed_energies + bs_res_full[point].nuclear_repulsion_energy\n",
    "    energies.append(energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxxUlEQVR4nO3deZQdVbn///enpySdqZN0ZjIBCRCmgA2CAyKTARlVFL6CgCjg7HJCf3i9Xq8oXJfXecpVVJRBHJAggSQgCDI3kIQMBEIIZE5n6sxJD8/vj1MNJ0130kmf03W6+/Naq1ZX7dpV9Zwjnie1d9XeigjMzMxyqSjtAMzMrOtxcjEzs5xzcjEzs5xzcjEzs5xzcjEzs5xzcjEzs5xzcrEuT9IvJf1HCte9V9Jl+3nsaElbJBXnOq5ck/RtSWslrWoet6SHJH0s7Rit48nvuVhnJmkJMBSoBxqA+cDNwJSIaEwxtH2SfI6PRcT9aceyLySNBhYCYyJiTQv7HwL+GBG/7ujYLF2+c7Gu4JyI6AuMAW4ArgV+k25IXYOkkr1UGQ2saymxWPfm5GJdRkTURsRU4EPAZZKOAJD0O0nfTtYrJf1D0kZJ6yU9Iqko2XetpOWSNktaKOnUpLyHpB9KWpEsP5TUo+m6ks6TNEvSJkkvS5qclL/eJCTpIEn/lLQuaUK6RVJFsu8PZH6k706alL4iaaykaPpxlzRC0tQk5kWSPp51/W9KukPSzUns8yRVtfY9Jef9rKTFSSzfy/oOLpf0qKQfSFoHfFNS/+TcNZJelfR1SUWSTgNmAiOSuH/XPO4Wrv1RSQskbZA0XdKY/fnf2gqfk4t1ORHxFLAMeGcLu7+Y7BtMpjnt/wNC0iHAp4Hjkrug9wBLkmOuA04AJgFHA8cDXweQdDyZZrgvAxXASVnHZRPwXWAEcBgwCvhmEu+lwGtk7sD6RMT/tHD87UncI4APAN+RdErW/nOTOhXAVOCnLZwj2wVAFXAscB7w0ax9bwUWk/l+rgd+AvQHDgTeBXwEuCJpwjsTWJHEffmeLijpPDLf9/vIfP+PALftJU7rpJxcrKtaAQxsobwOGE6mj6AuIh6JTMdjA9ADmCipNCKWRMTLyTEfBr4VEWsiogb4L+DSZN+VwE0RMTMiGiNieUS80PyiEbEoqbMzOcf/kvmh3itJo4C3A9dGxI6ImAX8msyPfJN/R8S0iGgA/kAmCe7JjRGxPiJeA34IXJy1b0VE/CQi6oFdwEXA1yJic0QsAb6f9fn3xTXAdyNiQXLu7wCTfPfSNTm5WFc1EljfQvn3gEXAjKRZ6KuQ+fEHPk/mbmKNpNsljUiOGQG8mnWOV5MyyNyBvMxeSBqanHO5pE3AH4HKNn6WEcD6iNjcLIaRWdursta3AT330l+ytNm5RrSyrxIo5c2fP/vabTUG+FHSJLmRzP8+2s9zWYFzcrEuR9JxZH6w/t18X/Kv7y9GxIFkmpK+0NS3EhG3RsQ7yPwIBnBjctiKpKzJ6KQMMj/EB7UhrO8k5zwyIvoBl5D5YX09tD0cuwIYKKlvsxiWt+G6rRnV7FwrsrazY1lL5m6v+effn2svBa6OiIqspVdEPLYf57IC5+RiXYakfpLOJtP38MeIeL6FOmdLOliSgFoyzWGNkg6RdErSUb8D2A40Pcp8G/B1SYMlVQLfIHPnAZmn0q6QdGrSyT1S0qEthNcX2ALUShpJpo8m22oyfRpvEhFLgceA70rqKekoMs1xf2ypfht9WdKApMntc8CfWrl2A3AHcL2kvkkT1hf289q/BL4m6XCA5EGBC/cvfCt0Ti7WFdwtaTOZfxlfR6Y/44pW6o4H7ifzQ/848POIeJBMf8sNZP6lvgoYAnwtOebbQDUwB3geeDYpa3p44ArgB2SS1b/Y/V/5Tf6LTOd5LXAP8Ldm+79LJoFtlPSlFo6/GBhL5g7jTuA/2/lOzF3AM8CsJJ49Pbr9GWArmU7+fwO3Ajft6wUj4k4yd4O3J02Dc8k8EGBdkF+iNOtmJAUwPulnMssL37mYmVnOObmYmVnOuVnMzMxyLpU7F0kDJc2U9FLyd0Ar9W6UNDdZPpRVfosyw3PMlXSTpNKk/GRJtcoMxTFL0jc66jOZmdkbUrlzkfQ/ZF4KuyF5iW1ARFzbrM57ybzUdiaZJ3keAk6NiE2SzgLuTareCjwcEb+QdDLwpYg4e1/iqaysjLFjx+7/BzIz64aeeeaZtRExuKV9exvxNF/OA05O1n9PJnFc26zORDJJox6olzQHmAzcERHTmipJego4oD3BjB07lurq6vacwsys25H0amv70urQHxoRK5P1VWQGyGtuNjBZUnny4tq72f2tYpLmsEuB+7KKT5Q0W5mJmg5vLQBJV0mqllRdU1PTrg9jZma7y9udi6T7gWEt7LoueyMiInnunmblM5JhPB4Dasi88NbQrNrPydzdPJJsP0tmQMItSdPZ38m8NPcmETEFmAJQVVXlpxrMzHIob8klIk5rbZ+k1ZKGR8RKScOBFicaiojryQz5jaRbgRezzvGfZIbtvjqr/qas9WmSfi6pMiLWtvsDmZlZm6XVLDYVaJpb/DIyQ1HsRlKxpEHJ+lHAUcCMZPtjZObbuDh7KltJw5Ixo5rm2SgC1uXxc5iZWQvS6tC/AbhD0pVkhu/+IIAys+ddExEfIzPM9yNJrtgEXJJ07kNmALxXgceT/X+LiG+RmUTpE5LqyQw8eFH4RR4zsw7nlyjJ9Ln4aTEzs30j6ZmIaHFKbQ//YmZmOZdWs5iZdZDGxmBXQ2Nmqc8sdQ2N1DUEDY1BXUMj9Y1BQ+MbZQ2NQUMEjY1vbDcGNEQQETRG0NAIja9vN61DRBDJdQOIZB/JepC9TovlvF4eNG9cyd7crW7Wnn1tkEm1/Sbl1qMJw/py9lEj9l5xHzm5mBWYuoZGarfXsXFbHbXb66jdvuv17S076tmyq55tOxvYurOeLTvr2bargW276tlR18iO+gZ21jWyo64hs9Q30tDopu9CJ+29Tr6cfdQIJxezziwi2LitjqUbtrF0/XaWb9xGzeadmWXLztfXN2yr2+N5ykqK6F1WTO8eJfQuK6F3j2LKy0oY2LuIHqXF9CotpmdpET1LiulRWkRZcTFlJUWUFoseJUWUFmeWkmJRWlxEcZEoLRYlRUWUFIniZCkqEsV6Y7u4SBQJiqQ3liKQMvWkzI+kyNSThMjU5/V9b5TDG/Wb1mmlvPlvr7Iqi92Pa6mOdTwnF7Mc21XfyMs1W1i4ajMvrNrMyzVbWLp+G8s2bGfLzvrd6vYoKWJIvx4M7tODcZW9OX7cQCr79GBg7zL69yqlf69SKsrfWO/bs4TSYneVWuFzcjFrh7qGRp5fXkv1kvU8v3wTC1dtYnHNVuqTpqjSYjFmUG9GDyznreMGMmpgOQcMKGfUwF4cMKCcfj1L/C9s65KcXMz2wbZd9Tz76kaeWrKep19Zz3NLN7CjLvMe78iKXhwyrC+nHTaUQ4b15dBh/RhX2ZuyEt9pWPfj5GK2Fxu37WLm/NXcN3cVjyxay676RooEhw3vx0XHjeb4cQOpGjuAIX17ph2qWcFwcjFrwZrNO5gxL5NQHl+8jobGYGRFLz781tG8a8Jgjh0zgH49S9MO06xgObmYJSKCJ19Zz+8eXcKM+atoDDiwsjdXn3Qgk48YxpEj+7t/xKyNnFys29tR18DUWSv47WNLWLByExXlpVx10kG879iRjB/SxwnFbD84uVi3tX7rLn7z78Xc+uRrbNhWx6HD+nLj+4/kvEkj6VlanHZ4Zp2ak4t1O7vqG7n58SX86IGX2LqzntMOG8oVbx/HCQcO9F2KWY44uVi3ERH884U1XH/PAhav3cq7JgzmP84+jIOH9E07NLMux8nFuoUXV2/mv/8xn0deWstBg3vz2yuO492HDEk7LLMuK7XkImkg8CdgLLAE+GBEbGih3o3Ae5PN/46IPyXlvwPeBdQm+y6PiFnJTJQ/As4CtiXlz+bvk1ghq2to5PszXuT/HllM77JivnH2RC49cYyHUDHLszTvXL4KPBARN0j6arJ9bXYFSe8FjgUmAT2AhyTdGxGbkipfjoi/NDvvmcD4ZHkr8Ivkr3UzKzZu5zO3Pcczr27gQ1WjuPbMQxnYuyztsMy6hTSTy3nAycn674GHaJZcgInAw8n0xvWS5gCTgTv2ct6bk+mNn5BUIWl4RKzMZfBW2B5cuIYv/GkWu+ob+fHFx3Du0bkfUtzMWpdm28DQrB/8VcDQFurMBiZLKpdUCbwbGJW1/3pJcyT9QFKPpGwksDSrzrKkzLqB+oZGbrzvBa747dMM7deTuz/zDicWsxTk9c5F0v3AsBZ2XZe9EREh6U0zGkXEDEnHAY8BNcDjQEOy+2tkklIZMIXMXc+39iG2q4CrAEaPHt3Ww6yArardwWdve46nlqzn4uNH8Z/nHO73VcxSktfkEhGntbZP0uqm5ipJw4E1rZzjeuD65JhbgReT8qa7np2Sfgt8Kdlezu53NwckZc3PO4VMUqKqqspT9XVyi9Zs4cO/foLNO+r54Ycmcf4xvlk1S1OazWJTgcuS9cuAu5pXkFQsaVCyfhRwFDAj2R6e/BVwPjA367wfUcYJQK37W7q2has2c9GUx2loDP76ibc5sZgVgDQ79G8A7pB0JfAq8EEASVXANRHxMaAUeCR5a3oTcEnSuQ9wi6TBZGY5nQVck5RPI/MY8iIyjyJf0SGfxlIxf8UmLvnNk5QUiVs/fgIHD+mTdkhmBijzUFX3VlVVFdXV1WmHYfvo+WW1XPKbJykvK+bWj5/AuMreaYdk1q1IeiYiqlra5zf0rVN67rUNfOSmp+jXs5TbrzqBUQPL0w7JzLL4NWXrdJ5esp5Lf/MUA8rL+NPVTixmhch3LtapLFy1mctveoqh/Xpy68dPYFh/Ty1sVoh852KdxsZtu/j4zdWU9yhxYjErcE4u1ik0NAafue05VtZu55eXvMWJxazAuVnMOoX/ue8FHnlpLTe870jeMmZA2uGY2V74zsUK3l2zlvOrhxdz6QljuOh4D9Vj1hk4uVhBm7u8lmv/Oofjxw7kP86emHY4ZtZGTi5WsNZt2cnVf3iGAeVl/OzDx1JW4v9czToL97lYQapvaOTTtz7H2i07+fM1JzK4b4+9H2RmBcPJxQrSrx5ezOOL1/H9C4/mqAMq0g7HzPaR2xms4Lxcs4UfPfASZx05jPe/5YC0wzGz/eDkYgWlsTH46l/n0Ku0mG+ee3ja4ZjZfnJysYJyy1Ov8fSSDVz33sMY0tcvSpp1Vk4uVjBW1m7nxntf4B0HV3Khm8PMOjUnFysIEcHX75xLQ2PwnQuOJJkgzsw6qVSSi6SBkmZKein52+J4HpJulDQ3WT6UVf6IpFnJskLS35PykyXVZu37Rgd9JGunu+es5IEX1vDFMyYwepCH0Dfr7NK6c/kq8EBEjAceSLZ3I+m9wLHAJOCtwJck9QOIiHdGxKSImAQ8Dvwt69BHmvZFxLfy+zEsF9Zv3cU3p87j6FEVXPH2cWmHY2Y5kFZyOQ/4fbL+e+D8FupMBB6OiPqI2ArMASZnV0iSzSnA3/MWqeXdf/9jPpu213Hj+4+kuMjNYWZdQVrJZWhErEzWVwFDW6gzG5gsqVxSJfBuYFSzOueTuQPalFV2oqTZku6V1OqzrJKuklQtqbqmpmb/P4m1y8Mv1nDnc8v55MkHceiwfmmHY2Y5krc39CXdDwxrYdd12RsREZKieaWImCHpOOAxoIZM81dDs2oXA7/O2n4WGBMRWySdReaOZnxL8UXEFGAKQFVV1Zuub/nX2Bh8Z9oCxgwq51OnHJx2OGaWQ3lLLhFxWmv7JK2WNDwiVkoaDqxp5RzXA9cnx9wKvJh1jkrgeOCCrPqbstanSfq5pMqIWNvuD2Q5d/ecFbywajM/umgSPUqK0w7HzHIorWaxqcBlyfplwF3NK0gqljQoWT8KOAqYkVXlA8A/ImJH1jHDlDzDKul4Mp9vXV4+gbXLrvpGvj/jRQ4b3o9zjhqRdjhmlmNpDVx5A3CHpCuBV4EPAkiqAq6JiI8BpcAjSa7YBFwSEfVZ57goOU+2DwCfkFQPbAcuigg3eRWgP1Uv5bX12/jt5cdR5E58sy5H/u3N9LlUV1enHUa3sX1XAyd970HGDirnjqtP9AuTZp2UpGcioqqlfX5D3zrcbx97hZrNO/nK5EOdWMy6KCcX61C12+r45UMvc8qhQzhu7MC0wzGzPHFysQ71q4dfZtOOer50xiFph2JmeeTkYh1mzaYd3PToK5x79AgmjvALk2ZdmZOLdZif/HMR9Q3BF06fkHYoZpZnTi7WIV5dt5XbnnqNDx03irGVvdMOx8zyzMnFOsSPH1hESbH47KktjsZjZl2Mk4vl3araHdw1azkXHTeaof08dbFZd+DkYnl38+NLaIjgo56rxazbcHKxvNq2q55bn3qN90wc5hkmzboRJxfLq78+u5yN2+q48p2+azHrTpxcLG8aG4Ob/v0KRx/Qn6oxA9IOx8w6kJOL5c0/X1jDK2u3cuU7D/QYYmbdjJOL5c2v/72YEf17cuYRLU1IamZdmZOL5cXc5bU8sXg9l799LKXF/s/MrLtJ7f/1ki6UNE9SYzJJWGv1JktaKGmRpK9mlY+T9GRS/idJZUl5j2R7UbJ/bL4+w8Ztu5g+bxU1m3fm6xKd1k3/foXysmI+dNzotEMxsxSk+U/KucD7gIdbqyCpGPgZcCYwEbhY0sRk943ADyLiYGADcGVSfiWwISn/QVIvL15dt42r//AMc5ZtzNclOqVVtTuYOnsFH6waRf9epWmHY2YpSC25RMSCiFi4l2rHA4siYnFE7AJuB85Tpnf4FOAvSb3fA+cn6+cl2yT7T1WeepMryjM/nBu21eXj9J2WX5o0s0JvDB8JLM3aXpaUDQI2RkR9s/Ldjkn21yb1dyPpKknVkqpramr2K7iK8jIg0zxmGdt21XPLk35p0qy7K8nnySXdD7T0qNB1EXFXPq+9NxExBZgCUFVVFftzjr49SigS1G73nUuTvz67nNrtdXzML02adWt5TS4RcVo7T7EcGJW1fUBStg6okFSS3J00lWcfs0xSCdA/qZ9zRUWif69SNvjOBYCI4LePvsLRoyp4i1+aNOvWCr1Z7GlgfPJkWBlwETA1IgJ4EPhAUu8yoOlOaGqyTbL/n0n9vBhQXsZG97kA8PSSDSyu2colbx3tlybNurk0H0W+QNIy4ETgHknTk/IRkqbB630mnwamAwuAOyJiXnKKa4EvSFpEpk/lN0n5b4BBSfkXgNcfX86H/uWlbhZL3FG9lD49SnjvUcPTDsXMUpbXZrE9iYg7gTtbKF8BnJW1PQ2Y1kK9xWSeJmtevgO4MKfB7kFFr1LWbnGz2OYdddwzZyXnHzOC8rLU/rMyswJR6M1iBa+ivMx9LsA9c1ayva6BC6tG7b2ymXV5Ti7tVFFeSq37XPhT9VLGD+nDMaMq0g7FzAqAk0s7VfQqY/POeuoaGtMOJTUvrd7Mc69t5INVo9yRb2aAk0u7Nb2lv6kbd+rfUb2UkiJxwbEj917ZzLoFJ5d26u5DwNQ1NPK3Z5dz6mFDqOzTI+1wzKxAOLm0U9MQMLXbu2en/gML1rBu6y4+6I58M8vi5NJOFcmov931Rco/Vy9lSN8evGvC4LRDMbMC4uTSTt25WWz1ph08uHAN73/LAZR4QjAzy+JfhHbqziMj//XZZTQGbhIzszdxcmmn7joyckTw5+plHD92IOMqe6cdjpkVGCeXdmoaGbm79bk8vWQDr6zdygeP812Lmb2Zk0sODOiGQ8A0DVJ51pEtTddjZt2dk0sOdLeRkbfsrOeeOSs55+jhHqTSzFrUpuQi6fuSDs93MJ1VRTdrFrt//mq21zVwwTEHpB2KmRWott65LACmSHpS0jWS+uczqM6mu42MfPfsFQzv35MqzzZpZq1oU3KJiF9HxNuBjwBjgTmSbpX07v25qKQLJc2T1Cipag/1JktaKGmRpK9mld+SlM+VdJOk0qT8ZEm1kmYlyzf2J7591Z1GRt64bRcPv1TD2UcNp6jIg1SaWcva3OciqRg4NFnWArPJzAR5+35cdy7wPuDhvVzvZ8CZwETgYkkTk923JHEcCfQCPpZ16CMRMSlZvrUfse2z7jQy8vR5q6hrCM45ekTaoZhZAWtTb6ykHwDnAA8A34mIp5JdN0pauK8XjYgFyXn3VO14YFEy4yRJEjsPmJ/MTtkU21NAqo3/2SMjD+rigzfePXslYwaVc+RIt4yaWevaeucyBzg6Iq7OSixN3jTVcI6MBJZmbS9Lyl6XNIddCtyXVXyipNmS7t3TQwiSrpJULam6pqamXYF2lyFgajbv5LGX13LOUSM8b4uZ7VFbnyOdDRzS7AelFng1ImpbOkDS/UBLL0FcFxF37VOUrfs58HBEPJJsPwuMiYgtks4C/g6Mb+nAiJgCTAGoqqqK9gTRXUZGvnfuShoDN4mZ2V61Nbn8HDiWzB2MgCOAeUB/SZ+IiBnND4iI09oZ23Ig+/XvA5IyACT9JzAYuDrrmpuy1qdJ+rmkyohY285Y9qi7jIx89+wVTBjah0OG9U07FDMrcG1tFlsBHBMRVRHxFuAYYDFwOvA/eYrtaWC8pHGSyoCLgKkAkj4GvAe4OCJe70WXNEzJ7ZWk48l8vnV5iu913aFZbMXG7Ty9ZAPnHOW7FjPbu7YmlwkRMa9pIyLmA4c2dbbvK0kXSFoGnAjcI2l6Uj5C0rTkGvXAp4HpZN6zuSMrhl8CQ4HHmz1y/AFgrqTZwI+BiyKiXU1ebdEdRka+Z85KAM52k5iZtUFbm8XmS/oF0PTY8YeSsh7APv9zPSLuBO5soXwFcFbW9jRgWgv1Wow7In4K/HRf42mv7jAy8t1zVnDkyP4eAdnM2qStdy6XAYuAzyfLYuByMollv16k7Eq6+sjIS9ZuZc6yWs45enjaoZhZJ7HXO5fkZcZpEfFu4PstVNmS86g6oa48MvI9z2eaxN7r/hYza6O93rlERAPQ6PHE9qwrj4x89+wVVI0ZwMiKXmmHYmadRFv7XLYAz0uaCWxtKoyIz+Ylqk6oolcpa7d0vTuXF1dv5oVVm/mvcz0otpm1XVuTy9+SxVpRUV7Gopqu10L4j9krKBKc6UnBzGwftCm5RMTvJfUCRkfEPo8l1h1UlJeycWvXahaLCO6es5ITDxrEkL490w7HzDqRtk4Wdg4wi2QML0mTJE3NY1ydTlccGXneik28snarX5w0s33W1keRv0lmgMqNABExCzgwLxF1UtkjI3cV985dSXGReM/hbhIzs33T1uRS18IAlV3nn+g50BWHgLlv7ireOm4gA3qXpR2KmXUybU0u8yT9P6BY0nhJPwEey2NcnU5XGxl50ZrNvFyz1XctZrZf2ppcPgMcDuwEbgM2kXlT3xJdbWTk6fNWA3DG4UNTjsTMOqO2Pi22DbguWawFTc1iXSe5rOLoURUM7+8XJ81s37V1muMJwJeAsdnHRMQp+Qmr82lqFusKQ8As37idOctquXbyoWmHYmadVFtfovwzmWHufw005C+czqsrjYw8Y94qAN7jJjEz209tTS71EfGLvEbSyXWlkZGnz1vFhKF9OHBwn7RDMbNOqq0d+ndL+qSk4ZIGNi3tubCkCyXNk9QoqWoP9SZLWihpkaSvZpX/TtIryWRhsyRNSsol6cdJ/TmSjm1PnPuioguMjLxuy06eemW9nxIzs3Zp653LZcnfL2eVBe17kXIu8D7gV61VSIb7/xmZ6ZSXAU9LmprMhAnw5Yj4S7PDzgTGJ8tbgV8kf/OuoguMjPzAgjU0Bk4uZtYubX1abFyuLxwRCwCSKe9bczywqGk6ZUm3A+cB8/dwzHnAzcn0xk9IqpA0PCJW5iby1nWFkZHvm7eKkRW9OHxEv7RDMbNObI/NYpK+krV+YbN938lXUFlGAkuztpclZU2uT5q+fpBMudyWYwCQdJWkaknVNTU1OQm2oryMjZ34JcotO+v590trmXzEsL0lfTOzPdpbn8tFWetfa7Zv8t5OLul+SXNbWM7b50jf7GvAocBxwEDg2n05OCKmRERVRFQNHjw4B+F0/pGRH3xhDbsaGt0kZmbttrdmMbWy3tL2m0TEafsc0e6WA6Oytg9Iyshq5top6bdk3sPZ4zH5lj0ycmlxW5+VKBzT562isk8ZbxkzIO1QzKyT29svYLSy3tJ2PjwNjJc0TlIZmTupqQCShid/BZxP5gEBkv0fSZ4aOwGo7Yj+FujcIyPvqGvgwRfWcPrEoRQXuUnMzNpnb3cuR0vaROYupVeyTrLdrtmjJF0A/AQYDNwjaVZEvEfSCODXEXFWRNRL+jQwHSgGboqIeckpbpE0OIllFnBNUj4NOAtYBGwDrmhPnPvi9SFgttcxqE+PvdQuLI+9vJatuxrcJGZmObHH5BIRxfm6cETcCdzZQvkKMsmhaXsamYTRvF6LQ88kT4l9KneRtl3TEDAbO+G7LtPnrqZvjxLedlBl2qGYWRfQ+ToGClhnHRm5vqGRmQtWc8phQygr8X8SZtZ+/iXJoc46MnL1qxtYv3WXm8TMLGecXHKos46MPH3eKnqUFPGuCbl5JNvMzMklhzrjyMgRwcz5q3nHwZX07tHW0YDMzPbMySWHOuPIyAtXb2bZhu2cPtHD65tZ7ji55FhmCJjOk1xmzluNBKccNiTtUMysC3FyybGK8tJO9SjyzAWrmTSqgiF92/XakpnZbpxccqyiEzWLrardwZxltW4SM7Occ3LJsc40MvL9C1YDcIaTi5nlmJNLjnWmkZFnzl/N2EHlHOTpjM0sx5xccix7ZORCtmVnPY+/vI7TJw713C1mlnNOLjnWWUZGfvjFGnY1NHL6RL+Vb2a55+SSY9kjIxeymfNXM6C8lGNHV6Qdipl1QU4uOdYZRkaua2jkny+s4ZRDh1LSCSc1M7PC51+WHOsMIyNXL9lA7fY6P4JsZnmTSnKRdKGkeZIaJVXtod5kSQslLZL01azyRyTNSpYVkv6elJ8sqTZr3zc64OPspjOMjDxz/mrKSop453jP3WJm+ZHWSIVzgfcBv2qtgqRi4GfA6cAy4GlJUyNifkS8M6veX4G7sg59JCLOzk/Ye1foIyNHBDMXrPJAlWaWV6ncuUTEgohYuJdqxwOLImJxROwCbgfOy64gqR9wCvD3vAS6Hwp9ZOQXV29h6XoPVGlm+VXIfS4jgaVZ28uSsmznAw9ExKasshMlzZZ0r6TD8xzjmxT6yMgz568C4NRDPVClmeVP3tpFJN0PtPQSxXURcVcL5fvjYuDXWdvPAmMiYouks8jc0YxvJb6rgKsARo8enaNwMgp5ZOSZ85OBKvt5oEozy5+8JZeIOK2dp1gOjMraPiApA0BSJZmmswuyrrkpa32apJ9LqoyItS3ENwWYAlBVVRXtjHU3hToy8upNO5i9rJYvv+eQtEMxsy6ukJvFngbGSxonqQy4CJiatf8DwD8iYkdTgaRhSsYykXQ8mc+3rgNjBgp3ZOSmgSrd32Jm+ZbWo8gXSFoGnAjcI2l6Uj5C0jSAiKgHPg1MBxYAd0TEvKzTXATc1uzUHwDmSpoN/Bi4KCJyelfSFoU6MvLM+asZM6ic8UM8UKWZ5Vcqz6JGxJ3AnS2UrwDOytqeBkxr5Rwnt1D2U+CnOQt0P2WaxQrrzmXzjjoeW7SOy942xgNVmlneFXKzWKdV0auMzTvqqS+gkZH/lQxUecbhHqjSzPLPySUPmt7SL6R3XabPW01lnzKOHT0g7VDMrBtwcsmDQhsZeWd9Aw++sIbTDhtKcZGbxMws/5xc8qDQRkZ+YvF6tuys54zD/ZSYmXUMJ5c8KLSRkafPW0XvsmLedpAHqjSzjuHkkgeFNDJyY2Mwc/5qTj5kCD1Li9MOx8y6CSeXPHi9WawA+lxmLdtIzeadbhIzsw7l5JIHfXuUUFos1mzesffKeTZ93ipKisTJh3igSjPrOE4ueVBUJMZV9ublNVtTjSMimDFvNSceNIj+ST+QmVlHcHLJk/FD+vLSms2pxvByzRZeWbvVL06aWYdzcsmT8UP78Nr6beyoa0gthunzkoEqD3N/i5l1LCeXPBk/pC8RsGjNltRimDFvFZNGVTCsv+duMbOO5eSSJxOGZkYeTiu5rKzdzuxltX5KzMxS4eSSJ2MG9aakSLy4Op1+l/vnZ5rEzpjo/hYz63hOLnlSVlLEuMrevJTSncv0eas5aHBvDvbcLWaWAieXPBo/tA8vpXDnUrutjicWr/NTYmaWmtSSi6QLJc2T1Cipag/1bpK0RtLcZuUDJc2U9FLyd0BSLkk/lrRI0hxJx+b7s7Rm/JC+qTwx9uDCNdQ3Bmd4OmMzS0mady5zgfcBD++l3u+AyS2UfxV4ICLGAw8k2wBnAuOT5SrgF7kIdn+MH9qHxsi8b9KRps9bxZC+PTj6gIoOva6ZWZPUkktELIiIhW2o9zCwvoVd5wG/T9Z/D5yfVX5zZDwBVEganoOQ99mEoX0BeGl1xyWXbbvqeWhhDWccPpQiz91iZinpzH0uQyNiZbK+CmhqAxoJLM2qtywp242kqyRVS6quqanJS4BjkyfGOvJN/fsXrGF7XQNnHzWiw65pZtZcXpOLpPslzW1hOS+X14mIAGIfj5kSEVURUTV48OBchvO6spIixlb25sUOvHOZOmsFw/r15PixAzvsmmZmzZXk8+QRcVoeT79a0vCIWJk0e61JypcDo7LqHZCUpWL8kD68sKpj7lxqt9XxrxfXcNmJY90kZmap6szNYlOBy5L1y4C7sso/kjw1dgJQm9V81uHGD+3Lq+u2dsgTY/fNW0ldQ3DuJDeJmVm60nwU+QJJy4ATgXskTU/KR0iallXvNuBx4BBJyyRdmey6AThd0kvAack2wDRgMbAI+D/gkx3ygVoxfkjmibHFNfkffn/q7BWMHVTOkSP75/1aZmZ7ktdmsT2JiDuBO1soXwGclbV9cSvHrwNObaE8gE/lLtL2ef2JsTWbmTiiX96us2bTDh57eR2fOWU8kpvEzCxdnblZrFMYW1lOcZHy/jjyP+asJALOPdpNYmaWPieXPOtRUszYQeV5H8By6uwVTBzez2OJmVlBcHLpAOOH9M3r0PuvrdvGrKUb3ZFvZgXDyaUDTBjahyV5fGLs7jkrADjHTWJmViCcXDrAwUP70hjwytr8PDE2ddYKqsYMYGRFr7yc38xsXzm5dICmWSnz0e/ywqpNLFy9mfPcJGZmBcTJpQOMq+xNcZHy0u8yddYKiovEWUemMjanmVmLnFw6QI+SYsbk4YmxiODuOSt4+8GVDOrTI6fnNjNrDyeXDjJ+SJ+cT3n83NKNLF2/3e+2mFnBcXLpIBOG9uXVddvYWZ+7J8amzlpBWUkR7zncM06aWWFxcukgBw/pQ0Nj5OyJsYbG4J7nV3LKIUPo27M0J+c0M8sVJ5cO0jTGWK7mdpk5fzU1m3dy/jFvmgfNzCx1Ti4dZFxlb4oEi3LQqR8R/PJfLzN6YDmnT3STmJkVHieXDtKztJixg3IzK+XTSzYwa+lGPn7SgRR7UjAzK0BOLh3o4CF9eGlN++9cfvWvlxnUu4wL33JADqIyM8u9VJKLpAslzZPUKKlqD/VukrRG0txm5d+T9IKkOZLulFSRlI+VtF3SrGT5ZZ4/yj6ZMLQvS9r5xNiLqzfzwAtruOxtY+lZWpzD6MzMcietO5e5wPuAh/dS73fA5BbKZwJHRMRRwIvA17L2vRwRk5LlmlwEmyvjh2aeGFuydtt+n2PKw4vpVVrMpSeMyWFkZma5lUpyiYgFEbGwDfUeBta3UD4jIuqTzSeATtE+NH5I0xNj+9c0trJ2O3fNWs6HjhvFgN5luQzNzCynukKfy0eBe7O2x0l6TtK/JL2ztYMkXSWpWlJ1TU1N/qMEDhyceWJsf9/U/+2jS2gMuPId43IcmZlZbpXk68SS7geGtbDruoi4K0fXuA6oB25JilYCoyNinaS3AH+XdHhEbGp+bERMAaYAVFVVRS7i2ZuepcWMGdSb2Us37vOxtdvruPXJ1zj7qOGMGlie++DMzHIob8klIk7L17kBJF0OnA2cGhGRXHMnsDNZf0bSy8AEoDqfseyLc44ewY8feInHFq3lbQdXtvm4W598jS0767nqpAPzGJ2ZWW50ymYxSZOBrwDnRsS2rPLBkoqT9QOB8cDidKJs2SdPPojRA8v5+t/ntvmpsZ31Ddz06Cu8c3wlh4/on+cIzczaL61HkS+QtAw4EbhH0vSkfISkaVn1bgMeBw6RtEzSlcmunwJ9gZnNHjk+CZgjaRbwF+CaiHjTAwFp6llazLfOO5zFa7fyq3+1Le/9/bnl1GzeydUnHZTn6MzMckNJi1K3VlVVFdXVHdty9qlbn2Xm/NXM+PxJjK3s3Wq9xsbgtB/8i16lxfzjM+9A8hv5ZlYYJD0TES2+q9gpm8W6gm+cPZGy4iL+4665tJbgI4KfPriIxTVbufpdBzmxmFmn4eSSkqH9evLFMybwyEtruef5lW/av6OugS/cMZv/nfki5xw9grOOaOnBOzOzwuTkkqJLTxjDESP78a2757N5R93r5Ws27eCiKU9w53PL+eLpE/jxRZMoKfb/VGbWefgXK0UlxUVcf/6R1GzZyfdnvAjA3OW1nPezR1m4ajO/vORYPnPqeDeHmVmnk7f3XKxtjh5VwaUnjOHmx5cwsHcZP39oEQPLy/jLJ070Y8dm1mn5zqUAfOk9hzCoTw/+d+aLHD6iP3d9+h1OLGbWqfnOpQD061nKTy4+hsdeXsen3n0QPUo8lL6ZdW5OLgXihAMHccKBg9IOw8wsJ9wsZmZmOefkYmZmOefkYmZmOefkYmZmOefkYmZmOefkYmZmOefkYmZmOefkYmZmOefJwgBJNcCraceRQ5XA2rSDKBD+Lt7g7+IN/i7e0J7vYkxEDG5ph5NLFySpurXZ4bobfxdv8HfxBn8Xb8jXd+FmMTMzyzknFzMzyzknl65pStoBFBB/F2/wd/EGfxdvyMt34T4XMzPLOd+5mJlZzjm5mJlZzjm5dCGSbpK0RtLctGNJm6RRkh6UNF/SPEmfSzumtEjqKekpSbOT7+K/0o4pTZKKJT0n6R9px5I2SUskPS9plqTqnJ7bfS5dh6STgC3AzRFxRNrxpEnScGB4RDwrqS/wDHB+RMxPObQOJ0lA74jYIqkU+DfwuYh4IuXQUiHpC0AV0C8izk47njRJWgJURUTOXyj1nUsXEhEPA+vTjqMQRMTKiHg2Wd8MLABGphtVOiJjS7JZmizd8l+Vkg4A3gv8Ou1YujonF+vyJI0FjgGeTDmU1CRNQbOANcDMiOiu38UPga8AjSnHUSgCmCHpGUlX5fLETi7WpUnqA/wV+HxEbEo7nrRERENETAIOAI6X1O2aTSWdDayJiGfSjqWAvCMijgXOBD6VNK3nhJOLdVlJ/8JfgVsi4m9px1MIImIj8CAwOeVQ0vB24Nykn+F24BRJf0w3pHRFxPLk7xrgTuD4XJ3bycW6pKQT+zfAgoj437TjSZOkwZIqkvVewOnAC6kGlYKI+FpEHBARY4GLgH9GxCUph5UaSb2Th12Q1Bs4A8jZk6ZOLl2IpNuAx4FDJC2TdGXaMaXo7cClZP51OitZzko7qJQMBx6UNAd4mkyfS7d/DNcYCvxb0mzgKeCeiLgvVyf3o8hmZpZzvnMxM7Occ3IxM7Occ3IxM7Occ3IxM7Occ3IxM7Occ3KxLknSljbU+byk8hxe83xJE3N4vsfaceyW5O8ISX/ZQ70KSZ/c3+uYtcbJxbqzzwP7lFwkFe9h9/lAzpJLRLwtB+dYEREf2EOVCsDJxXLOycW6NEknS3pI0l8kvSDpFmV8FhhB5uXCB5O6Z0h6XNKzkv6cjEvWNOfFjZKeBS6U9HFJTyfzo/xVUrmktwHnAt9LXtg8SNIkSU9ImiPpTkkDkvM9JOkHkqolLZB0nKS/SXpJ0rezYt+StX5tMu/GbEk3tPA5xyWxP9/sHGOb5veRdHgyr8usJKbxwA3AQUnZ9yT1kfRA8h08L+m8rPMskPR/yZwwM5K3/ZF0sKT7k9ielXRQUv7l5Huao24+h0y3FBFevHS5BdiS/D0ZqCUzYGMRmREM3pHsWwJUJuuVwMNk5j0BuBb4Rla9r2Sde1DW+reBzyTrvwM+kLVvDvCuZP1bwA+T9YeAG5P1zwEryLxF3wNY1nT+rM9wJvAYUJ5sD2zh804FPpKsfyrr2LHA3GT9J8CHk/UyoFf2/qS8hMw8J03fySJASb16YFKy7w7gkmT9SeCCZL0nmbvBM4ApybFFwD+Ak9L+78JLxy0lracdsy7jqYhYBpAMOz+WzIRZ2U4g06T1aGZYMsrIJKImf8paPyK5O6gA+gDTm19QUn+gIiL+lRT9HvhzVpWpyd/ngXkRsTI5bjEwCliXVfc04LcRsQ0gIlqas+ftwPuT9T8AN7ZQ53HgumROk79FxEvJZ90tdOA7yei4jWTmwBma7HslImYl688AY5OxqUZGxJ1JbDuSz3EGmQTzXFK/DzCeTAK3bsDJxbqDnVnrDbT8373IjLl1cSvn2Jq1/jsys1rOlnQ5mbuj/Y2psVl8ja3E1xZ7HMspIm6V9CSZybKmSboaWNys2oeBwcBbIqIuGUG4Z7OYIfM99trD5QR8NyJ+tQ/xWxfiPhfrzjYDfZP1J4C3SzoYXh8xdkIrx/UFViozpP+HWzpfRNQCGyS9M9l3KfAv9s9M4IqmJ9skDWyhzqNkRvqlWUyvk3QgsDgifgzcBRzF7t8BQH8yc57USXo3MGZPgUVmls9lks5PrtEjiXM68NGsfquRkoa05cNa1+DkYt3ZFOA+SQ9GRA1wOXBbMnrw48ChrRz3H2T6GR5l96Hrbwe+LOm5pFP7MjId/HOASWT6XfZZZEaqnQpUJ816X2qh2ufITPb0PK1P5/xBYG5yjiOAmyNiHZmmwLmSvgfcAlQl5/kIbRua/1Lgs8nnfAwYFhEzgFuBx5Nz/YXdk5h1cR4V2czMcs53LmZmlnNOLmZmlnNOLmZmlnNOLmZmlnNOLmZmlnNOLmZmlnNOLmZmlnP/P3LIz85jQYd6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "plt.plot(points, energies)\n",
    "plt.title(\"Dissociation profile\")\n",
    "plt.xlabel(\"Interatomic distance\")\n",
    "plt.ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_surface = EnergySurface1DSpline()\n",
    "\n",
    "xdata = res.points\n",
    "ydata = res.energies\n",
    "energy_surface.fit(xdata=xdata, ydata=ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.1576697027335805, -0.9127536865212236)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAELCAYAAAD3HtBMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsAElEQVR4nO3deXxU5dn/8c/FHqIJq0BdiguPa5UqoFZ/lkFbBUEFJS5Y0VpRFKsPtbjytNXaShQBreIO2roFlQJuyDKKuxBEES2KC4oiohYREYTk+v0xJxggy5DMmTOZ+b5fr3nNzMnJOde4zDf3fZ9z3+buiIiIbKtGURcgIiINkwJERETqRAEiIiJ1ogAREZE6UYCIiEidNIm6gHRq166dd+7cOfwTffQRrF4N++8f/rlEREJWWlr6pbu333J7TgVI586dmTdvXvgnuuACKCmBdJxLRCRkZra0qu3qwgpDixbw/fdRVyEiEioFSBjy8mDduqirEBEJlQIkDC1aQFkZbNgQdSUiIqFRgIQhLy/xrFaIiGQxBUgYWrRIPGscRESymAIkDC1bJp6/+y7aOkREQqQACUN+fuJZASIiWUwBEobttks8K0BEJIspQMKgFoiI5AAFSBgUICKSAxQgYagIkDVroq1DRCRECpAwqAUiIjlAARIGBYiI5AAFSBh0FZbINikuLiYej2+2LR6P06dPn622n3vuuZx77rkp37ehHqO4uJjIuHvOPA466CBPi/Jy98aN3a+4Ij3nk6wxatQonz179mbbhgwZ4kOGDNls2+zZs713795b7Vvd9kw/xpAhQ7xdu3ab9p89e7a3a9fOR48evdX2goICLywsTPm+DfUYW/4zDgMwz6v4To38Sz2dj7QFiLt7QYH7RRel73ySkaoKhJq+dKv6Ik3Zl09Bgcefecb9u+/8uccf993btPFb/vIX37tNG3+hpMT944/9xQce8INat/a7L7/ce7Ru7a9MnOi+aJG/evfdfkh+vh+23Xb+2u23u5eW+mu33eZHFhb6/Rdd5L8uLPR5//iH+0sv+bybb/behYX+4LBhfmxhoZeOHes+Z46Xjh3rx7Rs6X3y833+mDHu8bjPv/FGP66gwOePHu3zR4/24wsK/J7TT/cTCgp8/g03uM+a5fNvuMFPKCjwCaef7v0LCvz1G27w12+4wfsH2wYUFPjr11/vPmOGv3799T6goMAnDBq0aftW24qL3Z95xl8vLvYTCwp84qBBfmIt27dl3yiOkdRj5co6/3esAEl3gHTq5H722ek7n0SuqrAYPXq05+fn1/olH3/6ad+nTRt/5Z57vPSmm/yUggIv6d/fh+fn+3vnnecfnHGG/yMvz1896CB/sHlzX3HEEe7HHutf//zn/mqTJv5px47+duPG/t2OO7rvsouva9PGvzLz75s18+/Byxs1SvzvrkfuPp56qs7/bStAPM0Bssce7qeckr7zSdpU16qoqvWwU9u2PmHECD+poMAn9+3ro/LyfFm/fu7HHeer9t7bP2jUyL9v1iypL4D1TZr4CvCvCwvd997b/cAD3Q87zJfsuqtPA39rr73cTz3V/cwz3YcM8Ve6dfMbwecceqj7lVe6//nPPqNnT78U/KmjjnIfM8b9ppvcb7nFp/Tp42eDP9avn/vEie733ed+//3+cP/+XgT+0Iknuj/yiPujj/r9Awf6ceD/KipynzrV/fHH3Z94wu879VQ/Gvze005zf/pp9+nT3Z95xicMGuS9wCecfrr7zJnus2b5Pb/5jfcEv/s3v3GPxxOtkTFj/LjCQr/7jDO8X2Ghl44Z4/7cc146Zoz3Kyz0u844w/sGLZrSsWO9b2Gh3zV4cGLbuHHuzz/vpePG+bGFhX7n4MGJ1s+4cVVu8xde8NKbbvI+hYV+55lnep/CQi+96aZqt2/LvlEcI6nH11/X+b95BUi6A6RrV/d+/dJ3Pkm5pINi1iz/WevWvuDvf/f3hg71fzZv7h/usot/Vt1f/e3bu++/v/uvfuUL9tvPx4DP7NnT/ZZb3B980Bdcd50fXVjo4847z/dr3dqfmzbNZ8+Y4e3atfORI0dW2W2VzPZt2Tfdx9iyP19jIBoDybhHWgPksMPce/VK3/kk5ar7Unv+0Uf9zWuu8XF5eb5k1139S7PNAuLbli19DnjpAQe4X3ONL7riCu9bWOg3XHihd2zbNqkv14b0BZaKYwwZMiTpsaJtGczP9IsHUnGMUaNGedgUIOkOkF//2v3gg9N3PqmXmlobu7Vp4w+cdJLf3aKFr9l5501BsbFRI58HPvfnP3f/xz/c58zxOZMnbxYK2/LlWt0XaaZ/gTWUL0GpOwVIugOkf3/3ffdN3/mkXrb8y//Fhx/2y/Lz/asDD/SNQVfUuqZN3fv0cR81ykvHjfOd2rattfXQsmVLHz169Fbnqu5LV1+kkokUIOkOkNNPd+/cOX3nk3p7bto0v3C77fz9zp29LGhlrNllFx+Xl+d3DR7snYLup+q6tqrrhlEoSENXXYA0ieoGxqyXn6870TNQcXEx3bt3JxaLbdr26oQJNL3zTo5YuJAj1qzhvTVrePaXv6TlWWfR75JLKHniCWKxGLvF4xQVFTFgwABKSko2HSMWi1FSUsLcuXM3O27Fz7bcJpItFCBh2W47BUgG6t69O0VFRYkAaNKEr//3fzm4tJSypk1Z3qsXZ73yCj0uvJDxt93GgJdeUlCI1EABEpb8fFi7FsrLoZGmHMsUsViM6VddxXdHHw0bNrDRjPd/9zuW9+lD/yFDKJk8OREGvXpRVFTEKaecstXvKyhEEvTNFpaKGXnXro22jhxV1eR8Lz/wAIt/9jMOvPhiujZvzsXA7Zdeyu533slL771XbWtDRKqmAAmLpnSPVEVXVTweh3Xr+Oj00/n5oEHssWQJH5x1Fvs0b07ByJHcdNddxONxRowYUWW31IgRIyL6BCKZL5IAMbM2ZjbDzN4LnltXs98oM3sreJxcaftEM/vQzBYEj65pKz5ZCpBIVbQg/t6/P1/svDOd77+fVUcdxSv33svB06Zx36RJXH311ZSUlPwYNCKyTaJqgVwGzHL3LsCs4P1mzOxY4ECgK3AwcImZFVTa5Y/u3jV4LAi/5G2kNUGitWEDsRkzeHr1atZ/+SX3nXoqHWfM4MWPPlJXlUiKRDWIfjzQM3h9L/AscOkW++wDzHH3jcBGM3sTOAYoSVON9aMWSNpsdWnuZ5+x6phjaLVwIQ82b85Hv/89YydMYOegq2pLGhgXqZuoWiAd3H158PpzoEMV+7wBHGNmLc2sHRADdq7082vN7E0zG2Nmzas7kZkNMbN5ZjZv5cqVKfsAtVKApM1m4x3PPcf6/faj6cKF/C4vj5889RRXFherq0okBKEFiJnNrDR+UflxfOX9grscfcvfd/dngCeBl4AHgZeBsuDHlwN7Ad2BNmzdeql8nDvcvZu7d2vfvn1KPltSKgJkzZr0nTNHVXRDTT7uOMpiMT5ZvZrrTzqJQcENgJX3UVeVSOqE1oXl7kdV9zMzW2Fmndx9uZl1Ar6o5hjXAtcGv/MA8G6wvaL1st7MJgCXpLT4VFALJH3cib34IrE1a3gGmPuHP/DnUaO22k1dVSKpFVUX1lRgcPB6MDBlyx3MrLGZtQ1e7w/sDzwTvO8UPBtwAvBW+CVvIwVIepSVwdChMHIkDzdvzstXXMHYe+5RV5VIGkQ1iH4dUGJmZwNLgSIAM+sGnOfuvwOaAs8nMoLVwOnBgDrA/WbWHjBgAXBeestPgq7CCl95OZxzDkyYwLi8PPZ//HFO7tWLI4466sfpStTiEAlNJAHi7l8BR1axfR7wu+D1OhJXYlX1+71CLTAV1AJJuc2utiovh/POgwkTuH+PPdj/jjuSmrNKRFJHc2GFpVkzaNJEg+gptGkixIcfJvboo3DnnYzJy6Pr7bdrckORCChAwqQp3VOqomXxSt++xNau5ea8PLo+/jixXpnfIBXJRgqQMClAUi722WfE1q7lPmDlH/6g8BCJkCZTDJPWBEmt55+n/KyzeLFJEz68/HLG33abrrYSiZACJExqgaTOhx+yoW9f3i8vp/yRR/jT3/6mu8tFIqYACZMCJDXWr4eBAyn74Qe+mjiR/3d8YjID3V0uEi0FSJjy83UVVh1ttiDU8OFQWsp7V17JnM8+22w/rdkhEh0FSJjUAqmzikt2F40cCbfeyscDB9Jr3Di6d+8edWkiEtBVWGFSgNRZLBZj2pgx/PSMM/h4p504ePZsSiZN0r0dIhlELZAw6Sqsuisr45DbbqNxs2YctmwZ55x/vsJDJMMoQMKkFkjdjRsHL77IJU2bctbIkYwfP15XW4lkGAVImPLz4fvvEzPGSvL+8x/KLr+cp5s146QpU7R2uUiGUoCEqWJCxbVro62jISkrgzPP5IcmTdj+gQc23WmuS3ZFMo8G0cNUeUbe7bePtpaG4pZb4NVXyXvgAQ478cTNfqQJEkUyi1ogYdKU7tvm889h5Eg4+mg45ZSoqxGRWihAwqRFpWq12Q2DI0bAunW8MmgQxddfH21hIlIrBUiY1AKpVcUNg/PHjoV//pOPBg6k3/DhumFQpAHQGEiYKgJE05lUKxaLMemBB2jRuzerCgv55dNP64ZBkQZCARImtUCS0vO996CsjAHffMPgkSMVHiINhLqwwqQAqd233/LDVVfxUpMm7HfVVbphUKQBUYCESYPotfrwwgtp9t//0nzcOK6+5hrdMCjSgChAwqQWSM1WrGDHBx/kiyOO4KDzzwd0w6BIQ6IxkDBpEL1m11xDs7Iydrjjjs0264ZBkYZBLZAwNW2aeKgFsrUlS+D22+Gcc2DPPaOuRkTqQAESNs3IW7Vrr4UmTeD//i/qSkSkjhQgYVOAbO2DD+Cf/4Rzz4VOnaKuRkTqSAESNi0qBWwxZcl110Hjxrz4i19QXFwcbWEiUmcKkLDl52sQnR+nLHnpoYdg4kSW9e7NCRdcoClLRBowXYUVNnVhAT9envufPn3oUVZG3zlzKHn0UV1tJdKAKUDClp8PX34ZdRUZIbbnnhy+YQMTyss5btgwhYdIA6curLCpBbLJx8OH06isjDXDhmnKEpEsoAAJmwbRAZjzxBMUlpTwZc+e/O/NN2vKEpEskFSAmNlPzeyo4HWemWl91mSpBQLAhttvp9CdDqNGAZqyRCQb1DoGYmbnAEOANsDuwE7AbcCR4ZaWJXQVFpSVceTChXDYYdCjx6bNmrJEpGFLpgVyAXAYsBrA3d8DdgizqGyw6b6H/HxYvx7KyojH47l538PkyfDRRzB8eNSViEgKJRMg6939h4o3ZtYE8PBKyg4V9z0sWb4cgDlPPUVRUVFu3vdw442w225w/PFRVyIiKZRMgDxnZlcAeWb2K2ASMK2+JzazgWa2yMzKzaxbDfsdY2aLzWyJmV1WafuuZvZqsP1hM2tW35pSqaKP//Z//QuA8wcPpqSkJPe6bF55BV5+GS6+GBo3jroaEUmhZALkMmAlsBA4F3gSuCoF534LGADMqW4HM2sM3AL0BvYBTjWzfYIfjwLGuPsewH+Bs1NQU0rFYjEOOeooAH53yim5Fx4AN98MBQVw1llRVyIiKVZrgLh7OfAv4Cp3P8nd73T3endhufs77r64lt16AEvc/YOgG+0h4HgzM6AX8Eiw373ACfWtKdXi8TiPzZ4NwOP3358zl6xuGv/54gt45BEYPJj43Lm5Of4jksVqDRAzOw5YADwdvO9qZlNDrqvCjsAnld4vC7a1BVa5+8Yttm/FzIaY2Twzm7dy5cpQi60sHo9TVFTEH66+GoAbrroqZ+57qBj/ef/KK+GHH3j15z/P3fEfkSyWTBfWn0i0BFYBuPsCYNdkDm5mM83srSoeaRtNdfc73L2bu3dr3759uk7L3LlzKSkp4cBevQDo2rlzztz3EIvFKHnwQZrccw8f/vSn9B0xIjfHf0SyXDJzYW1w928SvUabJNWF5e5H1amqH30K7Fzp/U7Btq+AVmbWJGiFVGzPGCNGjEi8+PjjxPM33xA76aSc+RKN/fADlJdTtHQpQ0eOzJnPLZJLkmmBLDKz04DGZtbFzG4GXgq5rgpzgS7BFVfNgFOAqcEYTBw4KdhvMDAlTTVtm8LCxPOqVZGWkW5fXnMNX5ix7xVXaN4rkSyVTIBcCOwLrAceAL4BLq7vic2sv5ktAw4FnjCz6cH2n5jZkwBB62IYMB14Byhx90XBIS4FhpvZEhJjInfXt6ZQbL89mME330RdSdq8/OCDtHnlFdaedhp/uvZazXslkqVq7MIKLqN9wt1jwJWpPLG7TwYmV7H9M6BPpfdPkrh0eMv9PiAxNpPZGjVKXMaaQy0Qv+surFEjOv/tb8Dm816pK0ske9QYIO5eFtzoV+juufMndKq1apU7AVJWxi8WL4ajj4Zddtm0WfNeiWSfZAbR1wALzWwGsGlaWXf/fWhVZZtWrXKnC2vGDPj0Uxg7NupKRCRkyQTIY8GjMs2FtS0KC3OnBTJhArRtC/36RV2JiIQsmQBp5e7jKm8ws4tCqic7tWoFS5dGXUX4vvoK/v1vGDoUmjePuhoRCVkyV2ENrmLbmSmuI7vlSgvkgQfghx8075VIjqi2BWJmpwKnAbttMXXJ9sDXYReWVXJlEP2ee+DAA+GAA6KuRETSoKYurJeA5UA7YHSl7d8Cb4ZZVNZp1QpWr4by8sRlvdno9ddhwQL4xz+irkRE0qTaAHH3pcGNfuvc/bk01pR9CgvBHb799sc707PNxInQrBmcemrUlYhImtT457C7lwHlZpal33pp0qpV4jnLLuXdNG37xo3w0EPQrx/xN97QtO0iOUL3gaRDRYCsWrXZzXUNXcW07TMvuYQDvviChQccQFFRESUlJVGXJiJpUNf7QGRbZOmEihVTlCzu3Zv/adGCY8aNo2TSJN1xLpIjag0Qd783HYVktSztwgKI9ejB+vJyJm7YwNl//KPCQySHJLMiYRcze8TM3jazDyoe6Sgua2RpCwRg0d//TvMNG2gyeLCmbRfJMclcUzoBGA9sBGLAfSTWSJdkZWkLJB6Ps7y4mHU77MDZ99yjadtFckwyAZLn7rMAc/el7v5n4Nhwy8oyWdoCWRSP06usjBZnnw2NGm02bbuIZL9kBtHXm1kj4D0zG0Zi6djtwi0ryzRrBnl5WRcgw9q3T9wcOWjQpm2atl0kdyTTArkIaAn8HjgI+A1Vz48lNcnGKd0fegj23x/23TfqSkQkAslchVXRH7EG0Cx5dZVtEyp+8gm89BJce23UlYhIRGqaTHEaNaz74e7HhVJRtsq2FsikSYnnk0+Otg4RiUxNLZAbgmcD7gR+F345WaxVq8R6Gdni4YfhoINg992jrkREIlLTZIqbJlA0szWaULGeCgvh/fejriI1PvwQXnsNRo2KuhIRiVCyc4trCdv6yqYurIruq6KiaOsQkUjVNAbSptLbxmbWmkR3FgDurkWltkU2DaI//DD06AGdO0ddiYhEqKYxkFISLY+K0Jhf6WcO7BZWUVmpVavEcq/r1kGLFlFXU3dLlsD8+TB6dO37ikhWq2kMZNd0FpL1Kk/p3rFjlJXUT8VU7QMHRluHiEQuS9dXzUANeDqTTQtHATzyCBx6KPElS7RwlEiOU4CkSwOeULFi4aiX778fXn+dJcHCUd27d4+6NBGJUDJzYUkqVO7CamAqJkl8um9fDgVOfvhhSh59VHNeieQ4BUi6VHRhNcAWCCRCZI/WrSldu5Zjhw1TeIjItndhmdk7wWNYGAVlrQbcAgF4qaSEnT/9lP/27KmFo0QEqEOAuPvewOHAh6kvJ4s14EH0eDzOtN/+FoCjxo/XwlEiAiS3pO2FwU2Em7j7V+7+RHhlZaH8fGjcuEF2Yc2dO5cRXbrAPvvAXntp4SgRAZJrgXQA5ppZiZkdY2ZW62/I1swS3VgNsAUy4qyzaP3mm3DiiZu2xWIxRowYEWFVIhK1WgPE3a8CugB3A2eSWJnwb2amaVi3VWFhg2yB8O9/J1YeHDAg6kpEJIMkNQbi7g58Hjw2Aq2BR8xMd5JtiwbaAuGxx2C33eCAA6KuREQySDJjIBeZWSlQDLwI/Mzdh5JY3vbEGn+5+mMONLNFZlZuZt1q2O8YM1tsZkvM7LJK2yea2YdmtiB4dK1LHWnXECdUXL0aZs2C/v0T3XAiIoFk7gNpAwxw96WVN7p7uZn1reN53wIGALdXt4OZNQZuAX4FLCMxDjPV3d8Odvmjuz9Sx/NHo1UreO+9qKvYNk89BRs2wAknRF2JiGSYZAJkHGw1vfu37r7B3d+py0krfq+W8fgewBJ3/yDY9yHgeODtmn4pozXELqx//xvat4dDD426EhHJMMmMgcwHVgLvAu8Frz8ys/lmdlCIte0IfFLp/bJgW4VrzexNMxtjZs2rO4iZDTGzeWY2b+XKlWHVmpyGNoi+fj088QQcf3ziEmQRkUqSCZAZQB93b+fubYHewOPA+cCt1f2Smc00s7eqeByfgrovB/YCupPoYru0uh3d/Q537+bu3dq3b5+CU9dDq1bw7bewcWO0dSTr2WcT9ar7SkSqkEwX1iHufk7FG3d/xsxucPdza/rL392PqmdtnwI7V3q/U7ANd18ebFtvZhOAS+p5rvSouBt99Wpo06bmfTPB5MmJGyCPPDLqSkQkAyXTAlluZpea2U+DxwhgRTDIXR5ibXOBLma2q5k1A04BpgKYWafg2YATSAzKZ76GNKV7eTlMmQK9ezfsFRRFJDTJBMhpJP76/zcwmUSr4DSgMVBUl5OaWX8zWwYcCjxhZtOD7T8xsycB3H0jMAyYDrwDlLj7ouAQ95vZQmAh0A74a13qSLvWwYwwXzeA5eRfew0+/1zdVyJSrRq7sIJWxjh3H1TNLkvqclJ3n0wijLbc/hnQp9L7J4Enq9ivV13OG7kddkg8r1gRbR01KC4upnv37sSmT4cmTaBPH+LxeGI+LE1dIiKV1Bgg7l4WdFs1c/cf0lVU1urQIfGcwQFSsfrgRy1bkv/LXxJfsICioiJKKtZCFxEJJDOI/gHwoplNBb6r2OjuN4ZWVbZqAAESi8WYNno0+YMH88Tee3NmEB5aQEpEtpRMgLwfPBoB24dbTpbLz088MjhAAA754gsAzp8+naEjRyo8RKRKtQaIu/8FwMxauvva8EvKch06ZHyArLrvPpY1bszgK65g/PjxxGIxhYiIbCWZyRQPNbO3gf8E7w8ws2pvIJRaZHiAvDB5MtsvXEj+oEFcffXVWn1QRKqVzGW8Y4Gjga8A3P0N4IgQa8puHTtmdIB889BDNAZ2/f3vAbT6oIhUK5kxENz9ky0mPiwLp5wc0KEDPP981FVU69gNG2DHHeHAAzdtUxeWiFQlmRbIJ2b2C8DNrKmZXULixj6piw4d4KuvMnM+rHXrYPp0OO44rf0hIrVKJkDOAy4gMRPup0DX4L3URYcO4A5RzwxcldmzYe3axOy7IiK1SOYqrC+B6u5El21V+V6QTp2irWVLU6bA9ttDz55RVyIiDUCtAWJm7YFzgM6V93f334ZXVhbL1JsJy8th2jQ4+mhoXu0kyyIimyQziD4FeB6YiQbP668iQD7/PNo6tlRaCsuXq/tKRJKWTIC0dPdqF2ySbZSpLZCpU6FRo8T07SIiSUhmEP1xM+tT+26SlO22g7y8zAyQww+Htm2jrkREGohkAuQiEiGyzsxWm9m3ZrY67MKyllnm3Y2+dCm8+Wbi8l0RkSQlcxWWJlBMtUwLkGnTEs8KEBHZBsnMhWVmdrqZjQze72xmPcIvLYtlWoBMnQp77QVdukRdiYg0IMl0Yd1KYunZ04L3a4BbQqsoF2RAgBQXFycmSPzmG3j2WTjuOOLxOMXFxZHWJSINRzIBcrC7XwCsA3D3/wLNQq0q23XoAF9+CWXRXRVdsfLgWzfeCBs2MH+nnSgqKqJ79+6R1SQiDUsyl/FuCNZGd9h0Y2F5qFVluw4dEjfuffnlj5f1plnFLLvv9O7Nri1b0vsvf6Fk0iRNmigiSUumBXITMBnYwcyuBV4A/hZqVdkuQ+4FiR1+OH3NmLR2Leeef77CQ0S2Sa0B4u73AyOAvwPLgRPcfVLYhWW1DAmQ12+6ibx162g+cCDjx4/XolEisk2SXQ/kPwQrEkoKZECAxONx3r3qKvZv1oxTJ0yg42uvUVRURElJiVoiIpKUZLqwJNU6dkw8Rxggc197jTNataLxr38N+flaeVBEtpkCJAoFBYkZbyMMkBG9e5P3+eeb3TwYi8UYMWJEZDWJSMOiAEmz4uJi4s8+m+jGCmbkjeT+iylTEtOq9OuX3vOKSNZQgKRZxf0Xq4MJFePxeDT3X0ydCgcf/GN3mojINlKApFnFWMPLH3zA8jfeiGbg+tNPYd48rf0hIvWiAIlALBaj3b77wooVDB06NP1XPU2dmnhWgIhIPShAIhCPx5mzeDE7mHHbrbem//6LKVNgjz0SEyiKiNSRAiTNKsY8jj37bBq789idd1JUVJS+EPnmG5g9G044ITGILiJSRwqQNJs7dy4lJSX8z+GHA3B4ly7pvf/iySdhwwbo3z895xORrJXUneiSOpvus3jxxcTzxx8T69MnfeMgjz2WuPLqkEPScz4RyVpqgURljz0Sz++9l75zfv89PPVUYvC8kf7Vi0j96FskKjvsANtvn94AmTkTvvsOBgxI3zlFJGtFEiBmNtDMFplZuZl1q2G/e8zsCzN7a4vtbcxshpm9Fzy3Dr/qFDNLLCEbcoBsWnkQYPJkKCzk2WC7iEh9RNUCeQsYAMypZb+JwDFVbL8MmOXuXYBZwfuGJw0BUnHn+7MzZ8LUqXzevTsDBw3SyoMiUm+RBIi7v+Pui5PYbw7wdRU/Oh64N3h9L3BC6qpLoy5dYOlS+OGH0E5Rcef72BNPhK++4vJXX9WU7SKSEg11DKSDuy8PXn8ORLMubH116ZJY2vbDD0M9TSwW47I99+R7YDetPCgiKRJagJjZTDN7q4pHSufPcHcnWK+9mjqGmNk8M5u3cuXKVJ66/rp0STyH3I0VnzWLnUtL+eh//oeb7r5bKw+KSEqEdh+Iux8V1rGBFWbWyd2Xm1kn4Isa6rgDuAOgW7du1QZNJNIQIPF4nOtPPJEny8vZ8f/+j5Kf/EQrD4pISjTULqypwODg9WBgSoS11F3bttCqVagBMnfuXG4/8sjEAlb9+mnlQRFJGUv0AKX5pGb9gZuB9sAqYIG7H21mPwHucvc+wX4PAj2BdsAK4E/ufreZtQVKgF2ApUCRu1c12L6Zbt26+bx580L4RPXQowcUFsKMGeEcv7wcdt45cZ7Jk8M5h4hkNTMrdfetbrmIZCoTd58MbPVt5u6fAX0qvT+1mt//CjgytALTqUuXH6c1CcNLL8Fnn0FRUXjnEJGc1FC7sLJHly7w8cewbl04xy8pgRYtoG/fcI4vIjlLARK1Ll3AHT74IPXHLiuDRx6BPn0S06aIiKSQAiRqYV6J9eKLsHy5uq9EJBQKkKiFGSAlJZCXB8cem/pji0jOU4BErXXrxOW8qQ6QjRt/7L7abrvUHltEBAVIZghjUsVnnoEVK+D001N7XBGRgAIkE4QRIBMnJlo2ffrUuquISF0oQDJBly6wbBmsXZua4/33vzBlCpx2GjRrlppjiohsQQGSCSoG0t9/v16H2bR41MMPJ6aIP/NM4vG4Fo8SkVAoQDLBvvsmnufPr9dhKhaP+ubmm2G//YivWkVRUZEWjxKRUChAMsG++ybGK+o5zXosFmPaDTdQ+PbbTO/YkaKTT9asuyISGgVIJmjUCHr2TARIPSe3PGTxYsrMOHPmTIYOHarwEJHQKEAyRa9eiTmx6jOlycaNrLvrLmY3acI5I0cyfvx4LR4lIqFRgGSKipZCPb7w3/rrX2mxciUdR47k6quvpqSkhKKiIoWIiIRCAZIp9toLOnasV4C0mjiR7zt14mdXXAGgxaNEJFQKkAxRfP31rNhnH5g9e9M4yDZdgvvaa+y0dCl5l14KjRtv2hyLxRgxYkQYJYtIjlOAZIju3btz3Wuvweefw+LFxOPxbbsEd+xYKCiA3/421DpFRCooQDJELBajaPx4AKYNH05RUVHyl+AuWwaTJsHZZ2vdDxFJGwVIBjl00CBWFRSw7qmntu0S3FtuSax9fuGF4RYoIlKJAiSDxJ99lunr19OnZUtuu/XW5K6eWrUKbr8dTjgBdt017BJFRDZRgGSIijGP/S++mPy1a3n8uuuqvQR305xXAH/6E6xaxdyjj9acVyKSVgqQDDF37lxKSkrY++KLoXlzerzwQrWX4FbMefXa3XfDLbfwab9+9LnySs15JSJpZV7PqTMakm7duvm8efOiLqN2f/wjjB4NCxbA/vtXuUt89myaHX00XZs2pWteHnc88oimLRGRUJhZqbt323K7WiCZ6IoroFUruPTSaneJrVzJYRs3Mvz77zn1ggsUHiKSdgqQTNS6dSJEnn46cWPhlj74gPUXXMAbjRvT8corNeeViERCAZKphg3jm1atWD10aOIS3UDprbeyet99Wfv11/xw66385a9/1ZxXIhIJBUimatGCT887j4J332XN7rvDH//IuxdeyD4XXMD6pk1ZPGEC3YcMATTnlYhEQ4Pomay8nMXDh/P5rbdyWHk5TcrK+GaffSh89llo3z7q6kQkR2gQvSFq1Ig9x45l1mWXUVhWxh2//S2FpaUKDxHJCAqQDBePxxk/fjx/GDmSK6dOJf7yy1GXJCICKEAyWsXd6SUlJVogSkQyjgIkg1XcnV5xj4cGy0Ukk2gQXUREaqRBdBERSSkFiIiI1IkCRERE6iSSADGzgWa2yMzKzWyrfrVK+91jZl+Y2VtbbP+zmX1qZguCR5/wqxYRkcqiaoG8BQwA5tSy30TgmGp+NsbduwaPJ1NZnIiI1K5JFCd193cAzKy2/eaYWed01CQiItumIY+BDDOzN4NurtbV7WRmQ8xsnpnNW7lyZTrrExHJaqEFiJnNNLO3qngcn4LDjwd2B7oCy4HR1e3o7ne4ezd379Zec0iJiKRMaF1Y7n5UiMdeUfHazO4EHg/rXCIiUrUG2YVlZp0qve1PYlBeRETSKKrLePub2TLgUOAJM5sebP+JmT1Zab8HgZeBPc1smZmdHfyo2MwWmtmbQAz43zR/BBGRnKe5sEREpEaaC0tERFIqp1ogZrYSWBq8bQd8GWE5UdBnzg259plz7fNC+j/zT919q8tYcypAKjOzeVU1ybKZPnNuyLXPnGufFzLnM6sLS0RE6kQBIiIidZLLAXJH1AVEQJ85N+TaZ861zwsZ8plzdgxERETqJ5dbICIiUg8KEBERqZOcCxAzO8bMFpvZEjO7LOp60qG6lR2zlZntbGZxM3s7WPnyoqhrCpuZtTCz18zsjeAz/yXqmtLFzBqb2etmlhOTqprZR8FUTgvMLNKpNXJqDMTMGgPvAr8ClgFzgVPd/e1ICwuZmR0BrAHuc/f9oq4nbMFkm53cfb6ZbQ+UAidk879nS6zOlu/ua8ysKfACcJG7vxJxaaEzs+FAN6DA3ftGXU/YzOwjoJu7R37zZK61QHoAS9z9A3f/AXgISMX6JBnN3ecAX0ddR7q4+3J3nx+8/hZ4B9gx2qrC5QlrgrdNg0fW/3VoZjsBxwJ3RV1LLsq1ANkR+KTS+2Vk+RdLrguWRP458GrEpYQu6MpZAHwBzHD3rP/MwFhgBFAecR3p5MAzZlZqZkOiLCTXAkRyiJltBzwKXOzuq6OuJ2zuXubuXYGdgB5mltXdlWbWF/jC3UujriXNDnf3A4HewAVBF3Ukci1APgV2rvR+p2CbZJlgHOBR4H53fyzqetLJ3VcBceCYiEsJ22HAccGYwENALzP7V7Qlhc/dPw2evwAmk+iaj0SuBchcoIuZ7WpmzYBTgKkR1yQpFgwo3w284+43Rl1POphZezNrFbzOI3GhyH8iLSpk7n65u+/k7p1J/L88291Pj7isUJlZfnBhCGaWD/yaCFdkzakAcfeNwDBgOomB1RJ3XxRtVeGrYWXHbHUY8BsSf5EuCB59oi4qZJ2AeLBK51wSYyA5cVlrjukAvGBmbwCvAU+4+9NRFZNTl/GKiEjq5FQLREREUkcBIiIidaIAERGROlGAiIhInShARESkThQgIiJSJwoQERGpEwWISA4ws5vNbL6ZdY+6FskeChCRLBdMebEDcC6Q9etlSPooQESSZGZ/NrNLgtcv1bBfKzM7P32V1czdvyMx1cmzwE3RViPZRAEiUgfu/osaftwKyJgAMbO2QEvgW2BjxOVIFlGAiNTAzK40s3fN7AVgz0rb1wTP+Wb2RLAW+VtmdjJwHbB7MInj9cF+/w4WAFpUsQiQmXU2s3fM7M5g+zPBTLqY2Rlm9mZw3H9WOu/pwdrnC8zs9mCZ5tpcBdwALAL2TdE/GhGaRF2ASKYys4NITBPelcT/K/NJrK9e2THAZ+5+bPA7hSRWP9wvWNypwm/d/esgIOaa2aPB9i7Aqe5+jpmVACea2eskvvR/4e5fmlmb4Nh7AycDh7n7BjO7FRgE3FfDZ+gM/AIYDhxOIkCq7X4T2RYKEJHq/T9gsruvBTCzqtaOWQiMNrNRwOPu/ryZta5iv9+bWf/g9c4kguNz4EN3XxBsLwU6A62BSe7+JYC7V6xnfyRwEIkAAsgjsXxtTf4KXO3ubmbvoBaIpJACRKQe3P1dMzsQ6AP81cxmsUWLwMx6AkcBh7r7WjN7FmgR/Hh9pV3LSIRCdQy4190vT6Y2M+sKDAAON7NbgnMuTOZ3RZKhMRCR6s0BTjCzvGAVuH5b7mBmPwHWuvu/gOuBA0kMVm9fabdC4L9BeOwFHFLLeWcDA4PBbyq6sIBZwElmtkPFdjP7afB6lpntuMVxRgHHuXvnYNW+A1ALRFJILRCRarj7fDN7GHiDRFfR3Cp2+xlwvZmVAxuAoe7+lZm9aGZvAU+RGM84L+hCWgy8Ust5F5nZtcBzZlYGvA6c6e5vm9lVwDNm1ig43wVm9gmwB1DR1YWZ9QJauvvMSsddYWbbmVmbSt1iInWmFQlFGjgz24/EIP3wqGuR3KIAERGROtEYiIiI1IkCRERE6kQBIiIidaIAERGROlGAiIhInShARESkThQgIiJSJ/8fbA70/LuYEDoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xdata, ydata, \"kx\")\n",
    "x = np.arange(min(xdata) - 0.25, max(xdata) + 0.25, 0.05)\n",
    "plt.plot(x, energy_surface.eval(x), \"r-\")\n",
    "plt.xlabel(r\"distance, $\\AA$\")\n",
    "plt.ylabel(\"energy, Hartree\")\n",
    "dist = max(ydata) - min(ydata)\n",
    "plt.ylim(min(ydata) - 0.1 * dist, max(ydata) + 0.1 * dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation of the molecular Vibrational Energy levels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Born-Oppeheimer approximation removes internuclear vibrations from the molecular Hamiltonian and the energy computed from quantum mechanical ground-state energy calculations using this approximation contain only the electronic energy. Since even at absolute zero internuclear vibrations still occur, a correction is required to obtain the true zero-temperature energy of a molecule. This correction is called the zero-point vibrational energy (ZPE), which is computed by summing the contribution from internuclear vibrational modes. Therefore, the next step in computing thermodynamic observables is determining the vibrational energy levels. This can be done by constructing the Hessian matrix based on computed single point energies close to the equilibrium bond length. The eigenvalues of the Hessian matrix can then be used to determine the vibrational energy levels and the zero-point vibrational energy \t\n",
    "\\begin{equation}\n",
    "{\\rm ZPE} = \\frac{1}{2}\\, \\sum_i ^M \\nu_i \\, ,\n",
    "\\end{equation}\n",
    "with $\\nu_i$ being the vibrational frequencies, $M = 3N − 6$ or $M = 3N − 5$ for non-linear or linear molecules, respectively, and $N$ is number of the particles.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we fit a \"full\" energy surface using a 1D spline potential and use it to evaluate molecular vibrational energy levels.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAELCAYAAAD3HtBMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4nklEQVR4nO3de5yMdfvA8c9lndZinC2lthCJbLFSaHc8ehw7KdtBz48SJTo8klLpQCVbIj2ldNh9OpvCE1GhnRIl6xgdFJVSiAohp93r98c9M3axrDWH3Znr/XrNa/Z77733fc0ac+39/d7f6yuqijHGGHOsykQ6AGOMMaWTJRBjjDHFYgnEGGNMsVgCMcYYUyyWQIwxxhRL2UgHEE61atXSpKSk0J/oxx9h+3Y488zQn8sYY0JsyZIlW1S19sHbYyqBJCUlsXjx4tCfaNAg8HggHOcyxpgQE5F1h9tuXVihULEi/P13pKMwxpiQsgQSCvHxsHt3pKMwxpiQsgQSChUrQm4u7NsX6UiMMSZkLIGEQny882xXIcaYKGYJJBQqVnSebRzEGBPFLIGEQqVKzvPOnZGNwxhjQsgSSCgkJDjPlkCMMVHMEkgoVK7sPFsCMcZEMUsgoWBXIMaYGGAJJBQsgRhjYoAlkFDwJ5AdOyIbhzHGhJAlkFCwKxBjTAywBBIKlkCMMTHAEkgo2F1YxhyTjIwMvF5vgW1er5du3bodsv2GG27ghhtuCPq+pfUYGRkZRIyqxsyjVatWGhZ5eapxcap33x2e85moMWbMGM3Ozi6wbcCAATpgwIAC27Kzs7Vr166H7FvY9pJ+jAEDBmitWrUC+2dnZ2utWrV07Nixh2yvWrWqulyuoO9bWo9x8O84FIDFepjP1Ih/qIfzcTwJJDUzVTOXZaqq6t79ezU1M1VfWfGKqqru3LtTUzNT9c2Vb6qq6ta/t2pqvzidMrS7qqpu3rlZUzNTdfo301VVdcNfGzQ1M1Xf++49VVX9aetPmpqZqnPWzlFV1bV/rNXUzFT96IePVFX1m83faGpmqi74aYGqqq7ctFJTM1N10fpFqqq6bMMyTc1M1WUblqmq6qL1izQ1M1VXblqpqqoLflqgqZmp+s3mb1RV9aMfPtLUzFRd+8daVVWds3aOpmam6k9bf1JV1fe+e09TM1N1w18bVFV1+jfTNTUzVTfv3KyqqlO+mqKpmam69e+tqqr65so3NTUzVXfu3amqqq+seEVTM1N17/69qqqauSxTUzNTA7/LSYsn6T/++49A++lFT2uXV7sE2uM/G68Xvn5hoP3Ygse05+SegfboT0brFW9dEWiP/Gik9p7SO9AekT1C+/6vb6B915y7tP/0/oH27R/crje9e1Ogfet7t+qt790aaN/07k16+we3B9r9p/fXu+bcFWj3/V9fHZE9ItDuPaW3jvxoZKB9xVtX6OhPRquqkxA6/KeDPrbgscD3z51wrp7e9/TAf/wur3bRpxc9HfggLdevnA55bYiqOh8Scf3iNP7ceM3Ozta9+/dqy3EttfJ5lXXs2LFaM7GmthzXUt9c+aZmZ2drjfo1tOGohlqlTRXNzs7WzTs3a8txLTU+OV5dVavqlHff1NQXO2jGaw9qwxo19KEHhmil68vquJfuVf3pJ/VkjdPK15fVYXf31jbVq+vk50dr6tMp+vwzw7VtQoImJ8Xr2Y+cpos+fFkXPfustmlUWU+/5wRt26iyLv7Pf3TRey/o2Q810vanVtY3Bg/Wc5smaKtRDfWb917VJePHa5smFbTGdWV0+th7VL1enZhxo9bsF6ezHrtXl44dq23PrKRN7qyjXRIr69LHH9fpnlHaauSp2rVOZc285hptc3a8thp5qs57bJReWrWq3nhje63VL04/zXhYdc4cfWj0lVqrX5w+3/sq7Vm1qj74SLq2Gnmq9qxaVTN799az2lXUNiMbqc6ercsyMjS5Q0VtcUc9vaxqVV2WkaHjswbq+aNO18uqVtWs3r21hbuCdhzVXJdlZOhlVavq5YOT9cQry+qyjAzV2bN14MP/1AZXlNOs3r31sqpVtf9D/9C+T/8zsH/3W5ppUs9yuiwjQ5dlZGjjHuW1422NA+e79blL9eqH2gfOd+ol5fRfD50fiC+pZzntfkuzwP4XjWqtTbuUD5yv26izdOSL/6c6e7bq7Nl6xVOpOvql6wLtnhPa62OZ/Z325s1aXJZAjjOBuEa7tN87/VRVdeeeneoa7dIbZ9yoqk6CcI126S2zblFV1XVb16lruOgdQ1qoqpMAXKNdeu+H96qq6oqNK9Q12qWjPh6lqqoLf16ortEuHTN/jKqqen/wqmu0S8d/Nl5VVWd9O0tdo136bM6zqqo69aup6hrt0qxlWarqfIC7RrsCCSxrWZa6Rrt06ldTVVX12Zxn1TXapbO+naWqzge0a7RLvT94VVV1zPwx6hrt0oU/L1RV1VEfj1LXaJeu2LhCVVXv/fBedY12BRLQHbPvUNdol67buk5VVW+ZdYu6RrsCCebGGTeqa7RLd+5xEkq/d/qpa7Qr8Lv815R/afVHqwfaV7x1hdbKqBVoX/rmpVr3sbqBdvfXumv9sfUD7c6vdNYGTzQItDtmddSk8UmBdoeXOmijCY0C7bbPt9WmTzUNtFs910qbP9080E5+NlmTn00OtJs/3VxbPXfgvdL0qaba9vm2gXajCY20w0sdAu2k8UnaMatj4OqhwRMNtPMrnVVVdezYscrtaNsJzs9nZ2er3CHafFTzwF+PtTJqadqE87VZjRq68KWXtOqDlTQpvax6Lr1UhyQkaJX7ymv6kEb6n/h4/TjlLK10F9rn+kTV7t31uzYttPJd6HU94/WruDj9pnFddQ0XHXJhRf1dRFfULauuu9B7OqIKuqIO6roLHdXBaS+s77THtHPa3pOd9vg2TntWQ6f9bCunPbWp085q6bTfPMNpv3mG085q6bSnNnXaz7Zy2rMaOu3xbZy292SnPaad015Y32mP6uC0V9Rx2ve6nfY3NZ32HZ2c9rqqTvuWzk57c7zTvrG7095Z1mn3u8hpK87jX5eg1YcdaF9xOVrrjgPtS9PRukMPtLtfhdYfcqDduTfa4LYD7Y7/hybdeqDd4Vq00c0H2m37oU0HHWi36o82H3ignXyD8/C3mw909vG3mw5yjuFvN7rZOYe/nXSrE4O/3eA2J0Z/u/4Q5zUoqL73nhZXYQnExkCKKDkxmfYntQegXFw5khOTaXdSOwAqlatEcmIy5zU4DwBXBRfJWyvQdntVAGpWqklyYjJtTmgDQJ2EOiQnJtO6fmsA6lepT3JiMmfXOxuAk1wnkZyYTHJiMgCnVj+V5MRkWtRtAUDjmo1JTkymWe1mADSp1YTkxGSa1GoCQLPazUhOTKZxzcYAtKjbguTEZE6tfmrgtSQnJnOS6yQAzq53NsmJydSvUh+A1vVbk5yYTJ2EOgC0OaENyYnJ1KxUE4C2J7YlOTEZVwUXAOc1OI/kxGQqlXNqgLU7qR3JicmUiysHQPuT2gdeC0CHkzsEXivA+SefH/hdAKSenBr4XQGkJaXR9sS2Bdr+37W/3a5BuwJt/7+Vv93h5A4F2ueffH6B86WenFognrSktALx5m+7trlIIinQbtegHSfryaxdu5b09HQaVWhEWlIaXq+XcY88Qqv409g5bRn/u/BCcrp3p9O+Olz/I6ypXZukTp04e+Xv9HplHl/+8QfnXHcdKWt3cfea/fSaNo2xO3dy9k976Zq9hv779tFwyTLO2hxHp/VlYMMG6ldIoNnWitRa/Td5jRuT2KY9yVqHdqd04LtWrXjnr/2csqsq55zXCx54gG+Sz6XqRqBuMowbR/07HiS5SmN2N2jD9cDXrTuS7GpCcr+74bXXWJ3SifIbYVOzNHj7bRqPGEf93Nq8uBleS0+nyYjxJNdqTpMHJvDKVVfx7Gaol1eHxo88Cx98QIvh40jMq8ODWyHrmmtIvnMsyXVb8nnHS3ADv53RieS6Lan/4mSWjRvHzG2VSNS6DI2rwtJx42gz+BEaJTTklrJVePH//o/sPyvRKKEhPzw4mgtdLnY3PJ+4zXGsHv0YfPIJJ6ReQ9zmOF7u/S96uFyc2P5qGlVuRA+Xixf69GH55ngaVz0N5s9n6YQJrNpUkUQ5ge4uF0snTCD1kltoUr053V0uXujbl7WbKtK0ZkuWTphAd5eLhFpns39TOZZOmADz53PaGd3Yu7EcL/TtS3eXiyZNOtO+ebfA/tWrtGDnRmf/pRMm8MeG8tSKb3rgfG16cXpSauB82zeU5/SGHQPx7dxYjupVWgT2P73eOWzYUCFwvtPrtiat0/Uwfz7Mn895Tf5BWtcbAu22DVNJ63GT0z7nHILucFklWh9hGwNRVU1OVr3wwqPvZ0qsw41HHLa//sMPtUX16rp89Gj9buBAfaVCBf3hpJP01zJlAn8JFnjUrq165pmqF1ygy5s313Ggc9PSVJ9+WvWNN3T5o49qZ5dLn7zxRm1evbp+PGOGZs+Zo7Vq1dIRI0Ycts+8KNuPZd9wH+Pg/nwbA7ExkBL3CGsCaddOtWPH8J3PBF1hH2qfTJmiX4wapU/Gx+uaU07RLSIFEsRflSrpPNAlLVuqjhqlX959t/ZwufTxm2/WxJo1i/ThWpo+wIJxjAEDBoRkML+k3zwQjGOMGTNGQ80SyHEmkGMeRP93dZ1yUWNVtUH00jCI3m5Uu8B/WP8guv9qo8rlFfSiQU30xYoVdUeDBtr3YnSEG91fpowuBu3Sr7qOzOiuOm+ezps2TctfXV473tcx8CFa/l/l9Yb/3qCqzn/48n3K68WjLw58uHZ5tYve+uqtgQ/Ss8efrZMWTwrEV294PU291fn9+d97w98Yrl27dtVZc2YVeO9Nnz1dawypoQ94HlDVA++9zjd31gEDBhR472VnZ2vaxWnaclzLAu+9luNaapv0NpqdnV3gvTdgwAC9/KbLC7z3Jk2fpDWG1NBJ0514/e+9F2e8qF27dtUJ0yYUeO/1GNxD6w2vV+C913JcSx0+enjMvvdCdQOHqmrPyT0L3MBRXIUlEBsDCZW4OFuRsBSpm1iX9PT0wL33v65ezewLL2T04sVcsWcPzb5ZTe/cXBJatIDWrdjQuTNJ1avzzogReHft4ofadfDu30/P/v1p164dF1xwAR6PhxEjRtCkSRMaNWoEgNvtplWrVny7+ls8Hg9utxuA0047DY/HQ8OGDalevXqB2E477TT69u1bYFuzZs2YNWsWqampBbaf3+F8WrRoQYvmLQpsHzRoEM8991yBbW63m5f/+zLVqlUrsL1atWpkjMkIxOb33HPPcf/99xfYlpKSQosWLUhJSSmwvU2bNsyaNYtWrVoV2D709qGcdtpph5xv4MCBmFLocFklWh9h7cK65hrVpKTwnc8ct49nzNCbK1fWtUlJmuvrjtpx0kn6ZHy8vtCnj9bzdT8V1rVVWDdMOLoYjAklCrkCKRvpBBa1EhJsJnoJlJGRQUpKSoG/rj/PzKTc889z/sqVnL9jB9/t2MFHqalUuvZaLhw6FM/Mmbjdbk71eklPT6dnz54Frh7cbjcej4ecnJxD/mp3u92HbDMmWlgCCZXKlS2BlEApKSmkp6c7CaBsWf749785Z8kScsuVY0PHjly7cCFtbr6Zic8+S89PP7VEYcwRWAIJlYQE2LUL8vKgjA01lRRut5sP7r2XnZ07w7597Bdh7fXXs6FbNy4dMADPtGlOMujYkfT0dK688spDft4ShTEO+2QrorSsNLKWZwGwL3cfaVlpvPrFqwDs2reLtKw0Jq+aDMC23dtIK/sqU08Hdu1iy64tpGWlMWP1DAA27thIWlYa7695H4Cft/1MWlYac7+fC8D3f35PWlYaH//4MQCrt6wmLSuNT3/+FIBVv60iLSuNnF9yAFi+cTlpWWks37gcgJxfckjLSmPVb6sA+PTnT0nLSmP1ltUAfPzjx6RlpfH9n98DMPf7uaRlpfHztp8BeH/N+6RlpbFxx0YAZqyeQVpWGlt2bQFg6tdTSctKY9vubQBMXjWZtKw0du3bBcCrX7xKWlYa+3L3AZC1PIu0rLTA7/L5Jc/T6eVOgfYzOc/Q9bWugfaTC5/kojcuCrQf//RxLvNcFmg/Ov9Rrnz7wAf7qI9Hcc3UawLt+7z30fqh1oEB8eFzhzNgxgA+e/11VrdoweuzbuPlbnAb8Nydd/LUZQnc9eXYwNXGoJmDmLlvZuBqY8CMAQyfOzxw/GvfuZb7vPcF2tdMvYZRH48KtK98+0oenf9ooH2Z5zIe//TxQPuiNy7iyYVPBtpdX+vKMznPBNqdXu7E80ueD7SP+b2XlcbUr6cC2HsvAu+9a9+5NtD2v/f8hs4eyqCZgwLt296/jdvevy3QHjRzEENnDw20g/3eCzZLIKHiv+qwbqyIqFWr1oG7qnL385fXy1m9e9NozRq2nnEGb8SVoeqIEUx44QV+Xv8zrVu3Pmy31LBhwyL0CowpBQ43sh7qB1ADmAN853uuXsh+Y4BVvscV+bZnAT8Ay32P5KKcN6x3YWVmqoLq2rXhO6cpIDs7Wy9wuXRTrVqqoBs6ddL5kycf9g6qcMzmNaa0ooTNA7kL+FBVGwMf+toFiEh34GwgGTgHGCoiVfPtcoeqJvsey0Mf8jGyNUEia98+3HPm8P727ezZsoWXr7qKxDlzWPDjj4UOjBtjjk2kBtEvBtJ8X/8X+Ai486B9mgHzVHU/sF9EvgC6AJ4wxXh8bFXCsDnk1txff2Vrly5UW7mSNypU4MdbbmF8ZiYNvN7DdknZwLgxxROpK5C6qrrB9/VGoO5h9lkBdBGRSiJSC3ADDfJ9/2ER+UJExolIhcJOJCIDRGSxiCzevHlz0F7AUVkCCRv/rblerxc+/pg9zZtTbuVKro+Pp/5773FPRgYej6fATHNjzPELWQIRkbkisuowj4vz7+frX9ODf15VZwOzgE+BN4DPgFzft4cDTYEUnPGUg69e8h9nkqq2VtXWtWvXDsprKxJ/AtmxI3znjFH+bqhpF11ErtvNz9u389jll9PbNwEw/z7WVWVM8ISsC0tVOxX2PRHZJCL1VHWDiNQDfivkGA8DD/t+5nXgW992/9XLHhHJBIYe7ucjyq5AwkcV94IFuHfsYDaQc/vtPDBmzCG7WVeVMcEVqS6s6UAf39d9gHcO3kFE4kSkpu/rM4Ezgdm+dj3fswCX4NylVbJYAgmP3FwYOBBGjGByhQp8dvfdjH/pJeuqMiYMIpVAHgUuEJHvgE6+NiLSWkRe8O1TDvhERL4CJgHX+AbUAV4TkZXASqAW8FCoAz7myVwfXOVMJNy50yZzhWoyV14e9O/PqK+fo3WvOOrMmsX9Dz9M14yudH2u6yETCf0iPZnLJhJGwXvPJ9YnEkbkLixV/R34x2G2Lwau9329G+dOrMP9fMeQBhgMcXHOs12BBI3/bivK4dTKvfFGyMxkxaU1qO4+G3dH522RlJREWlraYWtWGWOCR5wx7NjQunVrXbx4cfhOWK4cDB0Ko0eH75xRzOurhuuZPBn3lCnwzDOMi48n+d13A8nDGBN8IrJEVVsfvN2KKYaSlXQPKv+dVAt79MC9axdPWfIwJqIsgYSSJZCgc//6K+5du3gZ2Hz77ZY8jIkgK6YYSrYmSHB98gl5117LgrJl+WH4cCY++6zdbWVMBFkCCSW7AgmeH35gX48erM3LI+/tt7n/kUdsdrkxEWYJJJQsgQTHnj3Qqxe5e/fye1YWHS52ihnY7HJjIssSSCglJFgpk2LKyMg4cGUxZAgsWcJ399zDvF9/LbCfrdlhTORYAimiYk3mqrfVJhJSvMlc/4v/H+np6Xw5YgSPL3+GrkNOoOOTT5KSklLqJ3PZRMKS/d6ziYRFZwkklCpWtC6sYqpVqxYzxo3j5Icf5k+XC+9vmwus42GMiTybSBhKN90Eb70F4SwjHy1ycyE1lb8XL+a0PXu4dsQIRo4cGemojIlJhU0ktCuQULJB9OJ78klYsICh5cpx7YgRTJw40e62MqaEsQQSSgkJ8Pffzl/Tpui++Ybc4cN5v3x5Ln/nHUaOHGm37BpTAlkCCSV/SfdduyIbR2mSmwt9+7K3bFmqvP56YKa53bJrTMljpUxCKf+aIFWqRDaW0uLpp+Hzz4l//XXaXXZZgW/ZglDGlCx2BRJKtqjUsdm4EUaMgM6d4corj76/MSaiLIGEUuXKzrMlkEIVmDA4bBjs3s3C3r3JeOyxyAZmjDkqSyChZFcgR5WSkkJ6ejpLx4+HV17hx169uHDIEGfhKGNMiWZjIKHkTyBWzqRQbrebt15/nYpdu7LV5SL1/ffxvPWWjXUYUwpYAgkluwIpkrTvvoPcXHpu20afESMseRhTSlgXVhEVqx7RHwsA2LJ9o9UjKqwe0V9/McEzlHOvEZrfey8TJ05k8KuDo7oekdXCKiHvPawW1vGyBBJKFSs6zzYPpFA/3HwzZXf9TZnTT2fkqFF4PB6yMrPYsmVLpEMzxhyF1cIKpa1boXp1eOIJ+Pe/w3fe0mLTJvaedBJb27alzscfBzZ7vV5ycnKsTLsxJURhtbBsDCSUbBD9yEaNonxuLnUmTSqw2SYMGlM6WBdWKJUr5zxsEP1Qa9bAc89B//7QpEmkozHGFIMlkFCziryH9/DDULYs3Hff0fc1xpRIlkBCzRLIob7/Hl55BW64AerVi3Q0xphisgQSapUrWwLhoJIljz4KcXEsOO88MjIyIhuYMabYLIGEWkKCDaJzoGTJp2++CVlZrO/alUsGDbKSJcaUYnYXVqhZFxZwYD2Pb7p1o01uLj3mzcMzZYrdbWVMKWYJJNQSEsAmxQHgbtKE9vv2kZmXx0WDB1vyMKaUsy6sULMrkICfhgyhTG4uOwYPtjXOjYkClkBCzQbRAZg3cyYuj4ctaWn8+6mnbI1zY6JAkRKIiJwsIp18X8eLSMytz1rsgnYJCWzJ/SvmC9rd7O2DS5W6Y8bw/JLnefjnhwNrnMdaQTsrpmjFFP2ivpiiiPQH3gae8206EfhfyCKKNtaFBXl5VN+5A9q1gzZtApvdbrfVuzKmFDtqMUURWQ60AT5X1bN821aqaovQhxdc4SymmJGRQUpKCu558+CBB2D/frzz5sVmkcC334ZevWDKFOjZM9LRGGOOUWHFFIvShbVHVffmO1BZIHZK+BaTf97Dmg0bAJj33nukp6fH5ryHJ56AU0+Fiy+OdCTGmCAqSgL5WETuBuJF5ALgLWDG8Z5YRHqJyJcikicih2S2fPt1EZHVIrJGRO7Kt/0UEfnct32yiJQ/3piCyT/v4blXnb7qm/r0wePxxN6tqwsXwmefwW23QVxcpKMxxgRRURLIXcBmYCVwAzALuDcI514F9ATmFbaDiMQBTwNdgWbAVSLSzPftMcA4VW0E/An0C0JMQeV2u2nbyRmwu/7KK2MveQA89RRUrQrXXnv0fY0xpcpRE4iq5gGvAveq6uWq+rwGYRUqVf1aVVcfZbc2wBpV/d7XjfYmcLGICNARZ3Af4L/AJccbU7B5vV6mZmcD8O5rr8XMLauBule//eaMf/Tpgzcnx+peGRNlinIX1kXAcuB9XztZRKaHOC6/E4Cf87XX+7bVBLaq6v6Dth9CRAaIyGIRWbx58+aQBpuf1+slPT2d20eOBODxe++NmXkP/vGftffcA3v38vlZZ8Xu+I8xUawoXVj341wJbAVQ1eXAKUU5uIjMFZFVh3mEbTRVVSepamtVbV27du1wnZacnBw8Hg9nd+wIQHJSUmDeQ7Rzu9143niDsi+9xA8nn0yPYcNic/zHmChXlFpY+1R1m9NrFFCkLixV7XT0vY7oF6BBvvaJvm2/A9VEpKzvKsS/vcQI3Kr700/O87ZtuC+/PGY+RN1790JeHunr1jFwxIiYed3GxJKiXIF8KSJXA3Ei0lhEngI+DXFcfjlAY98dV+WBK4HpvjEYL3C5b78+wDthiunYuFzO89atEQ0j3LaMGsVvIpxx991W98qYKFWUBHIzcAawB3gd2AbcdrwnFpFLRWQ9cC4wU0Q+8G2vLyKzAHxXF4OBD4CvAY+qfuk7xJ3AEBFZgzMm8uLxxhQSVaqACGzbFulIwuazN96gxsKF7Lr6au5/+GGre2VMlDpiF5bvNtqZquoG7gnmiVV1GjDtMNt/Bbrla8/CuXX44P2+xxmbKdnKlHFuY42hKxB94QWkTBmSHnkEODAnJicnx7qyjIkiR0wgqprrm+jnUtXY+RM62KpVi50EkpvLeatXQ+fOcNJJgc1ut9uShzFRpihdWDuAlSLyoohM8D9CHVhJc1wVUetUJq3OrNioiDpnDq/W/IU09zqriOpj1XitGq9fpN97wVaUu7Cm+h75WS2sY1G1KuxfH+kowiMz01kDpWbNSEdijAmxolTjvVVVnzzattIgnNV4C7j4Yli3DpYvD/+5w+n336F+fRg4EMaPj3Q0xpggOZ5qvH0Os63vcUcUS1yu2BgDef112LvX6l4ZEyMK7cISkauAq4FTDypdUgX4I9SBRZVYGUR/6SU4+2xo2TLSkRhjwuBIYyCfAhuAWsDYfNv/Ar4IZVBRp1o12L4d8vKc23qj0bJlThfdf/4T6UiMMWFSaAJR1XW+iX67VfXjMMYUfVwuUIW//jowMz3aZGVB+fJw1VWRjsQYEyZH/HNYVXOBPBGJ0k+9MKlWzXmOstnogbLt+/fDm2/ChRfiXbHCyrYbEyOKchuvfx7IHGCnf6Oq3hKyqKKNP4Fs3Vpgcl1p5y/bPnfoUFr+9hsrW7YkPT0dj8cT6dCMMWFQ3Hkg5lhEaUFFf4mS1V27clrFinR58kk8b71lM86NiRFHTSCq+t9wBBLVorQLC8Ddpg178vLI2rePfnfcYcnDmBhSlBUJG4vI2yLylYh873+EI7ioEaVXIABfjh5NhX37KNunj5VtNybGFOWe0kxgIrAfcAMv46yRbooqSq9AvF4vGzIy2F2nDv1eesnKthsTY4qSQOJV9UOcsifrVPUBoHtow4oyUXoF8qXXS8fcXCr26wdlyhQo226MiX5FGUTfIyJlgO9EZDDO0rGVQxtWlClfHuLjoy6BDK5d25kc2bt3YJuVbTcmdhTlCuRWoBJwC9AK+BeHr48V1Y67pPa/cpmxx5nAHzUltd9+Fc48k8m6ykpqWzl3wMq5l/T3XrAV5S4sf3/EDsCq5BVX2TjYufPo+5UWe/bA54vgvocjHYkxJkIKLecuIjM4wrofqnpRYd8rqSJWzh3g3HOd9dFnz47M+YPtiSfg9tthzRpo2DDS0RhjQqiwcu5HugLxX/cI8DxwfSgCixnVqjnrZUSLyZOhVStLHsbEsCMVUwwUUBSRHVZQ8Ti5XLB2baSjCI4ffoBFi2DMmEhHYoyJoKLWFrclbI9XtWrRMw/krbec5/T0yMZhjImoIy0oVSNfM05EquN0ZwGgqrao1LGIplUJJ0+GNm0gKSnSkRhjIuhIYyBLcK48/Eljab7vKXBqqIKKStWqOcu97t4NFStGOpriW7MGli6FsWOPvq8xJqodaQzklHAGEvXyl3RPTIxkJMfHX6q9V6/IxmGMibgoXV+1BCrF5UwCC0cBvP02nHsu3jVrbOEoY2KcJZBwKcUFFf0LR3322muwbBlrfAtHpaSkRDo0Y0wEFaUWlgmG/F1YpYy/SOL7PXpwLnDF5Ml4pkyxmlfGxDhLIOHi78IqhVcg4CSRRtWrs2TXLroPHmzJwxhz7F1YIvK17zE4FAFFrVJ8BQLwqcdDg19+4c+0NFs4yhgDFCOBqOrpQHvgh+CHE8VK8SC61+tlxnXXAdBp4kRbOMoYAxRtSdubfZMIA1T1d1WdGbqwolBCAsTFlcourJycHIY1bgzNmkHTprZwlDEGKNoVSF0gR0Q8ItJFROSoPxGFjntNhv+6mZFcCbZuLXVrMpzurs+lycvZcpmzbsLUr6fy4LoHueGWGwBbk8HWA7H1QPxK+nsv2I6aQFT1XqAx8CLQF2dlwkdExMqwHquESqXyCoSFC53n7raSsTHmgELXAzlkR5GWOAtKdQG8QFtgjqoOC114wRXR9UDAKX9erx68+27kYiiOrl3h22+dMiaxeQFqTEwrbD2QooyB3CoiS4AMYAHQQlUH4ixve9kRf7jwY/YSkS9FJE9EDgkq335dRGS1iKwRkbvybc8SkR9EZLnvkVycOMKuNBZU3L4dPvwQLr3UkocxpoCizAOpAfRU1XX5N6pqnoj0KOZ5VwE9gecK20FE4oCngQuA9TjjMNNV9SvfLneo6tvFPH9kVKsG330X6SiOzXvvwb59cMklkY7EGFPCFCWBPAmHlHf/S1X3qerXxTmp/+eOMh7fBlijqt/79n0TuBj46kg/VKJVq1b6rkD+9z+oXdtZktcYY/Ipyl1YS4HNwLfAd76vfxSRpSLSKoSxnQD8nK+93rfN72ER+UJExolIhcIOIiIDRGSxiCzevHlzqGItGperdA2i79kDM2fCxRc7tyAbY0w+RUkgc4BuqlpLVWsCXYF3gZuAZwr7IRGZKyKrDvO4OAhxDweaAik4XWx3Frajqk5S1daq2rp27dpBOPVxqFYN/voL9u+PbBxF9dFHTrzWfWWMOYyidGG1VdX+/oaqzhaRx1X1hiP95a+qnQr7XhH9AjTI1z7Rtw1V3eDbtkdEMoGhlAb+2ejbt0ONGkfetySYNs2ZAPmPf0Q6EmNMCVSUK5ANInKniJzsewwDNvkGufNCGFsO0FhEThGR8sCVwHQAEannexbgEpxB+ZKvNJV0z8uDd95xbuEtzSsoGmNCpigJ5Gqcv/7/B0zDuSq4GogD0otzUhG5VETWA+cCM0XkA9/2+iIyC0BV9wODgQ+ArwGPqn7pO8RrIrISWAnUAh4qThxhV91XEeaPUrCc/KJFsHGjdV8ZYwp1xC4s31XGk6rau5Bd1hTnpKo6DScZHbz9V6BbvvYsYNZh9utYnPNGXJ06zvOmTZGN4wgyMjJISUnB/cEHULYsdOuG1+t16mENKzVzRo0xYXDEBKKqub5uq/KqujdcQUWtunWd5xKcQPyrD/5YqRIJqal4ly8nPT0dj38tdGOM8SnKIPr3wAIRmQ7s9G9U1SdCFlW0KgUJxO12M2PsWBL69GHm6afT15c8bAEpY8zBipJA1voeZYAqoQ0nyiUkOI8SnEAA2v72GwA3ffABA0eMsORhjDmsoyYQVX0QQEQqqequ0IcU5erWLfEJZOvLL7M+Lo4+d9/NxIkTcbvdlkSMMYcoSjHFc0XkK+AbX7uliBQ6gdAcRQlPIPOnTaPKypUk9O7NyJEjbfVBY0yhinIb73igM/A7gKquAM4PYUwlUtAW9UlMZOPW9SV2UZ9tb77JnEZwbetVbNyxEbfbzW0Tb+P6+dfboj4+tqCULSjlV9ree8FWpDXRVfXngzblhiCW2FC3LmzZEukoCtV93z6oWROqHBjuOrPFmTRo0OAIP2WMiUVHXVBKRN4GngD+A5wD3Aq0VtUrj/iDJVDEF5QCeOABGDkS9u515lmUJLt3O8mjTx94xnopjTGOYi8oBdwIDMKphPsLkOxrm+KoWxdUIdKVgQ8nOxt27XKq7xpjzFEU5S6sLUBhM9HNsco/F6RevcjGcrB33nG6rtLSIh2JMaYUOGoCEZHaQH8gKf/+qnpd6MKKYiV1MmFeHsyYAZ07Q4VCiywbY0xAUTrh3wE+AeZig+fHz59ANm6MbBwHW7IENmyw7itjTJEVJYFUUtVCF2wyx6ikXoFMnw5lyjjl240xpgiKMoj+roh0O/pupkgqV4b4+JKZQNq3d+7CMsaYIihKArkVJ4nsFpHtIvKXiGwPdWBRS6TkzUZftw6++AIuuujo+xpjjE9R7sKyAorBVtISyAxnlrIlEGPMsShKLSwRkWtEZISv3UBE2oQ+tChW0hLI9OnQtCk0bhzpSIwxpUhRurCewVl69mpfewfwdMgiigUlIIFkZGQ4BRK3bYOPPoKLLsLr9ZKRkRHRuIwxpUdREsg5qjoI2A2gqn8C5UMaVbTz18PKjdxd0f6VB1c98QTs28fSE08kPT2dlJSUiMVkjCldinIb7z7f2ugKgYmFeSGNKtrVretM3Nuy5cBtvWHmdrvxeDx83bUrp1SqRNcHH8Tz1lu27ocxpsiKcgUyAZgG1BGRh4H5wCMhjSralZC5IO727ekhwlu7dnHDTTdZ8jDGHJOjJhBVfQ0YBowGNgCXqOpboQ4sqpWQBLJswgTid++mQq9eTJw40RaNMsYckyLVE1fVb/CtSGiCoAQkEK/Xy7f33suZ5ctzVWYmiYsWkZ6ejsfjsSsRY0yRFGlBKRPkVeGqCGl94f1fnFXfIrEq3PWf9KNDg8rE/fOfzN30GQ+ue5D/vPwfcnJybFU4W5HQViT0ibb3XrBZAomEKlWcGelbt0YshAY1a1Jxy5YCkwfPO+88hg0bFrGYjDGly1FXJIwmJWFFwoyMDFJSUnD37QupqfDyy3i9XnJycsL74T1qFNx/P/z6KyQmhu+8xphS53hWJDRB5J9/sd1XUNHr9UZm/sX06XDOOZY8jDHFZgkkzPzzLz77/ns2rFgRmYHrX36BxYtt7Q9jzHGxBBIBbrebWmecAZs2MXDgwPDf9TR9uvNsCcQYcxwsgUSA1+tl3urV1BHh2WeeCf/8i3fegUaNnAKKxhhTTJZAwsw/5tG9Xz/iVJn6/POkp6eHL4ls2wbZ2XDJJc6dYMYYU0yWQMIsJycHj8fDae3bA9C+cWM8Hg85OTnhCWDWLNi3Dy69NDznM8ZErSLNRDfBE7hVd8EC5/mnn3B36xa+cZCpU507r9q2Dc/5jDFRy65AIqVRI+f5u+/Cd86//4b33nMGz8vYP70x5vjYp0ik1KnjzEgPZwKZOxd27oSePcN3TmNM1IpIAhGRXiLypYjkicghsxvz7feSiPwmIqsO2l5DROaIyHe+5+qhjzrIRJwlZEOcQAIrDwJMmwYuFx/5thtjzPGI1BXIKqAnMO8o+2UBXQ6z/S7gQ1VtDHzoa5c+YUgg/pnvH82dC9OnszElhV69e9vKg8aY4xaRBKKqX6vq6iLsNw/44zDfuhj4r+/r/wKXBC+6MGrcGNatg717Q3YK/8z38ZddBr//zvDPP7eS7caYoCitYyB1VXWD7+uNQGTWhT1ejRs7S9v+8ENIT+N2u7mrSRP+Bk61lQeNMUESsgQiInNFZNVhHkGtn6FOOeFCSwqLyAARWSwiizdv3hzMUx+/xo2d5xB3Y3k//JAGS5bw42mnMeHFF23lQWNMUIRsHoiqdjr6XsW2SUTqqeoGEakH/HaEOCYBk8Ap5x7CmI5dGBKI1+vlscsuY1ZeHifcdx+e+vVt5UFjTFCU1i6s6UAf39d9gHdCfcKQrAr3Zw5Uq8bPa5eFbFW4nJwc+l+YTNq1wup2TXG73dwz6R6un3+9rQrnYysS2oqEftH+3gu2SN3Ge6mIrAfOBWaKyAe+7fVFZFa+/d4APgOaiMh6Eenn+9ajwAUi8h3Qydcuffy38oZwDGTY0KHU/fJLqFkDKlcG4KyzzqJBgwYhO6cxJjbYioSR1ru3U9bkxx9Dc/z586FDB3j9dbjqqtCcwxgT1WxFwpKqcWP46SfYvTs0x/d4oGJF6NEjNMc3xsQsSyCR1rgxqML33wf/2Lm58Pbb0K2bUzbFGGOCyBJIpIXyTqwFC2DDBkhPD/6xjTExzxJIpIUygXg8EB8P3bsH/9jGmJhnCSTSqleHmjWDn0D27z/QfeW7+8oYY4LJEkhJEIqiirNnw6ZNcM01R9/XGGOKwRJISRCKBJKV5VzZdOsW3OMaY4yPJZCSoHFjWL8edu0KzvH+/BPeeQeuvhrKlw/OMY0x5iCWQEoC/0D62rXHdZjA4lGTJzsl4vv2xev12uJRxpiQsARSEpxxhvO8dOlxHca/eNS2p56C5s3xbt1Kenq6LR5ljAkJSyAlwRlnOOMVx1lm3e12M+Pxx3F99RUfJCaSfsUVVnXXGBMylkBKgjJlIC3NSSDHWZus7erV5IrQd+5cBg4caMnDGBMylkBKio4dnZpYx1PSZP9+dr/wAtlly9J/xAgmTpxoi0cZY0LGEkhJ4b9SOI4P/FUPPUTFzZtJHDGCkSNH4vF4SE9PtyRijAkJSyAlRdOmkJh4XAmkWlYWf9erR4u77wacMRGPx0NOTk6wojTGmABLICVExmOPsalZM8jODoyDHNMtuIsWceK6dcTfeSfExQU2u91uhg0bFoqQjTExzhJICZGSksKjixbBxo2wejVer/fYbsEdPx6qVoXrrgtpnMYY42cJpIhCvS71g+se5JSxzlrJzw+/kX++8U/umXQPbrf76OtSr/iAtIpvsLz/hVCliq1LXcLXpbY10e2952dropugafaPf7C1alX2fPwxJ9Q/gbPOOqtoPzh5Mii2ZK0xJqxsTfQSxOv18lvXrvSIi+OU+Hgmv/XW0edxbN0Kp57q3MU1ZUpY4jTGxBZbE72E8495nHnbbSTs2sW7jz5a6C24gZpXAPffD1u3ktO5s9W8MsaElSWQEiInJwePx8Ppt90GFSrQZv78Qm/B9de8WvTii/D00/xy4YV0u+ceq3lljAkr68Iqie64A8aOheXL4cwzD7uLNzub8p07k1yuHMnx8Ux6+20rW2KMCQnrwipN7r4bqlWDO+8sdBf35s2027+fIX//zVWDBlnyMMaEnSWQkqh6dSeJvP++M7HwYN9/z55Bg1gRF0fiPfdYzStjTERYAimpBg9mW7VqbB84EPLyApuXPPMM2884g11//MHeZ57hwYcesppXxpiIsARSUlWsyC833kjVb79lR8OGcMcdfHvzzTQbNIg95cqxOjOTlAHOBCWreWWMiQQbRC/J8vJYPWQIG595hnZ5eZTNzWVbs2a4PvoIateOdHTGmBhhg+ilUZkyNBk/ng/vugtXbi6TrrsO15IlljyMMSWCJZASzuv1MnHiRG4fMYJ7pk/H+9lnkQ7JGGMASyAlmn92usfjsQWijDEljiWQEsw/O90/x8MGy40xJYkNohtjjDkiG0Q3xhgTVJZAjDHGFIslEGOMMcUSkQQiIr1E5EsRyRORQ/rV8u33koj8JiKrDtr+gIj8IiLLfY9uoY/aGGNMfpG6AlkF9ATmHWW/LKBLId8bp6rJvsesYAZnjDHm6MpG4qSq+jWAiBxtv3kikhSOmIwxxhyb0jwGMlhEvvB1c1UvbCcRGSAii0Vk8ebNm8MZnzHGRLWQJRARmSsiqw7zuDgIh58INASSgQ3A2MJ2VNVJqtpaVVvXthpSxhgTNCHrwlLVTiE89ib/1yLyPPBuqM5ljDHm8EplF5aI1MvXvBRnUN4YY0wYReo23ktFZD1wLjBTRD7wba8vIrPy7fcG8BnQRETWi0g/37cyRGSliHwBuIF/h/klGGNMzLNaWMYYY47IamEZY4wJqpi6AhGRzcA6X7MWsCWC4USCvebYEGuvOdZeL4T/NZ+sqofcxhpTCSQ/EVl8uEuyaGavOTbE2muOtdcLJec1WxeWMcaYYrEEYowxplhiOYFMinQAEWCvOTbE2muOtdcLJeQ1x+wYiDHGmOMTy1cgxhhjjoMlEGOMMcUScwlERLqIyGoRWSMid0U6nnAobGXHaCUiDUTEKyJf+Va+vDXSMYWaiFQUkUUissL3mh+MdEzhIiJxIrJMRGKiqKqI/Ogr5bRcRCJaWiOmxkBEJA74FrgAWA/kAFep6lcRDSzEROR8YAfwsqo2j3Q8oeYrtllPVZeKSBVgCXBJNP87i7M6W4Kq7hCRcsB84FZVXRjh0EJORIYArYGqqtoj0vGEmoj8CLRW1YhPnoy1K5A2wBpV/V5V9wJvAsFYn6REU9V5wB+RjiNcVHWDqi71ff0X8DVwQmSjCi117PA1y/keUf/XoYicCHQHXoh0LLEo1hLICcDP+drrifIPlljnWxL5LODzCIcScr6unOXAb8AcVY361wyMB4YBeRGOI5wUmC0iS0RkQCQDibUEYmKIiFQGpgC3qer2SMcTaqqaq6rJwIlAGxGJ6u5KEekB/KaqSyIdS5i1V9Wzga7AIF8XdUTEWgL5BWiQr32ib5uJMr5xgCnAa6o6NdLxhJOqbgW8QJcIhxJq7YCLfGMCbwIdReTVyIYUeqr6i+/5N2AaTtd8RMRaAskBGovIKSJSHrgSmB7hmEyQ+QaUXwS+VtUnIh1POIhIbRGp5vs6HudGkW8iGlSIqepwVT1RVZNw/i9nq+o1EQ4rpEQkwXdjCCKSAPyTCK7IGlMJRFX3A4OBD3AGVj2q+mVkowq9I6zsGK3aAf/C+Yt0ue/RLdJBhVg9wOtbpTMHZwwkJm5rjTF1gfkisgJYBMxU1fcjFUxM3cZrjDEmeGLqCsQYY0zwWAIxxhhTLJZAjDHGFIslEGOMMcViCcQYY0yxWAIxxhhTLJZAjDHGFIslEGNigIg8JSJLRSQl0rGY6GEJxJgo5yt5UQe4AYj69TJM+FgCMaaIROQBERnq+/rTI+xXTURuCl9kR6aqO3FKnXwETIhsNCaaWAIxphhU9bwjfLsaUGISiIjUBCoBfwH7IxyOiSKWQIw5AhG5R0S+FZH5QJN823f4nhNEZKZvLfJVInIF8CjQ0FfE8THffv/zLQD0pX8RIBFJEpGvReR53/bZvkq6iMj/icgXvuO+ku+81/jWPl8uIs/5lmk+mnuBx4EvgTOC9KsxhrKRDsCYkkpEWuGUCU/G+b+yFGd99fy6AL+qanffz7hwVj9s7lvcye86Vf3DlyByRGSKb3tj4CpV7S8iHuAyEVmG86F/nqpuEZEavmOfDlwBtFPVfSLyDNAbePkIryEJOA8YArTHSSCFdr8ZcywsgRhTuA7ANFXdBSAih1s7ZiUwVkTGAO+q6iciUv0w+90iIpf6vm6Akzg2Aj+o6nLf9iVAElAdeEtVtwCoqn89+38ArXASEEA8zvK1R/IQMFJVVUS+xq5ATBBZAjHmOKjqtyJyNtANeEhEPuSgKwIRSQM6Aeeq6i4R+Qio6Pv2nny75uIkhcII8F9VHV6U2EQkGegJtBeRp33nXFmUnzWmKGwMxJjCzQMuEZF43ypwFx68g4jUB3ap6qvAY8DZOIPVVfLt5gL+9CWPpkDbo5w3G+jlG/zG34UFfAhcLiJ1/NtF5GTf1x+KyAkHHWcMcJGqJvlW7WuJXYGYILIrEGMKoapLRWQysAKnqyjnMLu1AB4TkTxgHzBQVX8XkQUisgp4D2c840ZfF9JqYOFRzvuliDwMfCwiucAyoK+qfiUi9wKzRaSM73yDRORnoBHg7+pCRDoClVR1br7jbhKRyiJSI1+3mDHFZisSGlPKiUhznEH6IZGOxcQWSyDGGGOKxcZAjDHGFIslEGOMMcViCcQYY0yxWAIxxhhTLJZAjDHGFIslEGOMMcViCcQYY0yx/D+5i+fzkvLTpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vibrational_structure = VibrationalStructure1DFD(mol, energy_surface)\n",
    "\n",
    "plt.plot(xdata, ydata, \"kx\")\n",
    "x = np.arange(min(xdata) - 0.25, max(xdata) + 0.25, 0.05)\n",
    "plt.plot(x, energy_surface.eval(x), \"r-\")\n",
    "plt.xlabel(r\"distance, $\\AA$\")\n",
    "plt.ylabel(\"energy, Hartree\")\n",
    "dist = max(ydata) - min(ydata)\n",
    "plt.ylim(min(ydata) - 0.1 * dist, max(ydata) + 0.1 * dist)\n",
    "for N in range(15):\n",
    "    on = np.ones(x.shape)\n",
    "    on *= energy_surface.eval(\n",
    "        energy_surface.get_equilibrium_geometry()\n",
    "    ) + vibrational_structure.vibrational_energy_level(N)\n",
    "    plt.plot(x, on, \"g:\")\n",
    "on = np.ones(x.shape)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a partition function for the calculation of heat capacity\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The partition function for a molecule is the product of contributions from translational, rotational, vibrational, electronic, and nuclear degrees of freedom. Having the vibrational frequencies, now we can obtain the vibrational partition function  $q_{\\rm vibration}$ to compute the whole molecular partition function \n",
    "\\begin{equation}\n",
    "q_{\\rm vibration} = \\prod_{i=1} ^M \\frac{\\exp\\,(-\\Theta_{\\nu_i}/2T)}{1-\\exp\\,(-\\Theta_{\\nu_i}/2T} \\, . \n",
    "\\end{equation} \n",
    "Here $\\Theta_{\\nu_i}= h\\nu_i/k_B$, $T$ is the temperature and $k_B$ is the Boltzmann constant.  \n",
    "\n",
    "The single-point energy calculations and the resulting partition function can be used to calculate the (constant volume or constant pressure) heat capacity of the molecules. The constant volume heat capacity, for example, is given by \n",
    "\n",
    "\\begin{equation}\n",
    "C_v = \\left.\\frac{\\partial U}{\\partial T}\\right|_{N,V}\\, ,\n",
    "\\qquad\n",
    "{\\rm with} \\quad\n",
    "U=k_B T^2 \\left.\\frac{\\partial {\\rm ln} Q}{\\partial T}\\right|_{N,V} .\n",
    "\\end{equation}\n",
    "\n",
    "$U$ is the internal energy, $V$ is the volume and $Q$ is the partition function. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we illustrate the simplest usage of the partition function, namely creating a Thermodynamics object to compute properties like the constant pressure heat capacity defined above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjoUlEQVR4nO3deZRb9X338fdXo9kXjz0eL3jBC2BwXDBhIDZNApTlodmTQ1JIF0JISNskJWmfpw2nfRra05zTPE/SbE9LShKSNIEkDQ2FkoYlLoEQltQGA8YL2NgYYxuPx+vMeKSR9H3+uFdjYY9saeZqpJE+r3N0pLvo3q+uNPc7v/v73d/P3B0REZHRxModgIiIVC4lCRERyUtJQkRE8lKSEBGRvJQkREQkLyUJERHJa8KThJndZmZ7zGxdzrz3m9nzZpYxs56JjklEREZXjpLEd4Arj5m3Dngf8MiERyMiInnFJ3qH7v6ImS04Zt4GADOb6HBEROQEJjxJjJeZ3QDcANDa2nremWeeWeaIREQmlzVr1ux19+5C1p10ScLdbwVuBejp6fHVq1eXOSIRkcnFzF4udF21bhIRkbyUJEREJK9yNIH9AfA4sMTMdpjZ9Wb2XjPbAawEfmpm9090XCIicrxytG66Js+iuyY0EBEROSldbhIRkbyUJEREJC8lCRERyUtJQkRE8lKSEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8lCRERCQvJQkREclLSUJERPJSkhARkbyUJEREJC8lCRERyUtJQkRE8lKSEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8lCRERCQvJQkREclLSUJERPJSkhARkbyUJEREJC8lCRERyWvCk4SZ3WZme8xsXc68aWb2oJm9GD5Pnei4RETkeOUoSXwHuPKYeZ8BVrn76cCqcFpERMpswpOEuz8C7Dtm9ruB74avvwu8ZyJjEhGR0VVKncRMd98Vvt4NzCxnMCIiEqiUJDHC3R3wfMvN7AYzW21mq3t7eycwMhGR2lMpSeI1M5sNED7vybeiu9/q7j3u3tPd3T1hAYqI1KJKSRL3ANeGr68F7i5jLCIiEipHE9gfAI8DS8xsh5ldD/w9cLmZvQhcFk6LiEiZxSd6h+5+TZ5Fl05oICIiclKVcrlJREQqkJKEiIjkpSQhIiJ5KUmIiEheShIiIpKXkoSIiOSlJCEiInkpSYiISF5KEiIikpeShIiI5KUkISIieSlJiIhIXkoSIiKSl5KEiIjkpSQhIiJ5KUmIiEheShIiIpKXkoSIiOSlJCEiInkpSYiISF5KEiIikpeShIiI5KUkISIieSlJiIhIXkoSIiKSl5KEiIjkVXSSMLO/KEUgIiJSeeInW8HM/jV3ElgOfL5UAU2UZ3ccYM+hBJctnVnuUEREKlYhJYlD7v6B8PF+4OelCsbMbjSzdWb2vJl9qlT7AfjcTzfw0e+t5mfP7SrlbkREJrVCksTnjpn+y1IEYmbLgI8CFwDnAO8ws9NKsa/hdIZndhzAgE/9aC3rXj1Yit2IiEx6J00S7r4VwMymh9P7ShTLWcCT7j7o7ingYeB9pdjRhl2HGBrOcPO73kDMjDvX7CjFbkREJr1iKq5vK1kUgXXAW8ysy8xagLcB845dycxuMLPVZra6t7d3TDta8/J+AC5fOpPfmDuFta8cGHvUIiJVrJgkYSWLAnD3DQQV4g8A9wFrgfQo693q7j3u3tPd3T2mfa15eT+nTGli9pRmls/rZP3OQyRSx+1KRKTmFZMkvGRRZHfg/i13P8/d3wrsB14oxX6eenk/bzx1KgDL53WSTGfYsOtwKXYlIjKpVUxJAsDMZoTP8wnqI+6Ieh87Dxxh58EhenKSBMDa7fuj3pWIyKR30vskctxUsiiO+jcz6wKGgY+7+4God/D8zkMAnBMmh9lTmpjR3qh6CRGRURSUJMxsnruvy7PsHe5+bxTBuPtbotjOifQeTgAwe0ozAGbG8nmdPLNDzWBFRI5V6OWmB81swbEzzezDwFcijajE+vqDJDGttWFk3jnzOtm6d4CDg8PlCktEpCIVmiT+FHjAzE7PzjCzm4BPAxeVIrBS2dufoKMpTkP86Ec/a3Y7AJt7VXktIpKroMtN7v6fZpYAfmZm7wE+QnBn9FvdfVLV+O4dSDK9vfF18xZ3twGwZc8A5506rRxhiYhUpIJbN7n7KuA64BfAIuC3JluCgOBy0/TW1yeJuVNbaKiLsaW3v0xRiYhUpkIrrg8T3CdhQCNwKbDHzAxwd+8oXYjR2tuf5PQZba+bVxczFkxvYUvvQJmiEhGpTIVebmovdSATpa8/wYpFx19SWtzdxqbdqpMQEclVUyPTpdIZ9g8OM72t8bhli7vbeHnfIMlUpgyRiYhUpppKEvsGkwB0jZYkZrSSzjjb9+mSk4hIVk0lib2HgyQxPeceiaxsC6fNe5QkRESyaipJ9A0EN9Id2wQWYOH0VgBe2qsWTiIiWUUnCTP7i1IEMhH6+sPLTaOUJNqb6pnZ0cgWlSREREactHWTmf1r7iSwnGDch0lnb9glx2h1EgCLprfpXgkRkRyFNIE95O4fyU6Y2S0ljKek9vYnaaiL0dE0+sde2N3Kfz63a4KjEhGpXIVcbvocHB3jGvjL0oVTWn39CbraGgjuATzewq5WDgwOcyBsBSUiUutOmiTcfWv48rZwel9JIyqhvWGSyGdBWHm9da/qJUREoMJGpiu1voEkXa2j10cALJzeAsC2PiUJERGosDGuS62vPznq3dZZ86a1EDPYqj6cRESAGitJ7BtIMrWlPu/yxngdc6Y2s7VvcAKjEhGpXMUkiYkY47pkUukMR4bTtDflTxIAC7pa2aY6CRERoLjxJEYd43qyGEikAWjL0/w1a9H0IEm4T/qrayIi41bseBIjszg6vsSkGE/icCIYv7q98cQfecH0Vg4nUuztT9I9SvcdIiK1pGbGk+hPpICTlySyzWC39Q0oSYhIzRtL303nmNknwsfZpQiqFPqHwiRxkpLEIt0rISIyoqgkYWY3ArcDM8LH7Wb2yVIEFrXDBZYk5nQ2E4+ZKq9FRCjwclOO64E3ufsAgJl9Hngc+FrUgUUtW5I4WZ1EvC7G/GktKkmIiFD85SYD0jnTaSbJ/ROF1klAUC+hJCEiUnxJ4tvAk2Z2Vzj9HuBbkUZUIoXWSUAwANHjW/pw97ydAYqI1IKiShLu/g/Ah4F94eM6d/9yVMGY2afN7HkzW2dmPzCzpqi2na2TaG0orCRxZDjNa4cSUe1eRGRSKrYkgbuvAdZEHYiZzQH+BFjq7kfCwY6uBr4Txfb7h1K0NcaJxU5eMljYdXQo01lTIstTIiKTTrGtm3rM7C4ze8rMnjWz58zs2QjjiQPNZhYHWoCdUW24PzFMewH1ERAMPgSwba/6cBKR2lZsSeJ24H8BzwGZKANx91fN7AvAduAI8IC7P3DsemZ2A3ADwPz58wvefn8iVVB9BMDsjiYa4zF1GS4iNa/Y1k297n6Pu29195ezjygCMbOpwLuBhcApQKuZ/d6x67n7re7e4+493d3dBW//8FCqoJZNALGYcWqXmsGKiBRbkvismX0TWAWM1Oq6+08iiOUyYKu79wKY2U+AC4HvR7DtokoSELRw2qJxJUSkxhWbJK4DzgTqOXq5yYEoksR2YIWZtRBcbroUWB3BdoGg4np2EZXQC6a38tDGXtIZp66Aym4RkWpUbJI4392XlCIQd3/SzO4EngJSwNPArVFtv9iSxOLuNpLpDK/sGxzp9E9EpNYUWyfxmJktLUkkgLt/1t3PdPdl7v777h7ZjQpBE9gTDziU67QZbQBs3tMfVQgiIpNOsUliBbDWzDaVqAlsSWQyTn+y8IprOJokXlSSEJEaVuzlpitLEkWJDQ6ncT955365OprqmdnRqJKEiNS0opJEVM1dJ9pIv01FlCQgKE1s3nO4FCGJiEwKRQ86NBn1h0OXFlNxDXBadxtbejXetYjUrppIEofHUZLoT6TYfWioFGGJiFS8mkgS2bEkiqmTADhtRjC094uvqV5CRGpTQUnCzA6b2aFRHofN7FCpgxyv8dRJgJrBikjtKuis6e7tpQ6klEbGty6yJDG9rYEpzfVqBisiNavo8STCjvhOB0b6uHD3R6IMKmpHx7cu/GY6ADNjyax2Nu2u+MKSiEhJFDuexEeAR4D7gb8Jn2+OPqxoZeskWhvrin7v0tkdbNx9mExGLZxEpPYUW3F9I3A+8LK7XwKcCxyIOqio9SdSNNfXEa8rvp7+zFntDCbTbN+nAYhEpPYUe9YccvchADNrdPeNQEk6/ItSfyI1plIEwFmzOwDYsEuXnESk9hSbJHaYWSfw78CDZnY3UPF3YQ8mUrQWWWmdtWRWOzFTkhCR2lRstxzvDV/ebGYPAVOA+yKPKmIDyTQtDWNLEk31dSyc3sr6XeqeQ0Rqz9jOnIC7PxxlIKU0mEzR2jC2y00QXHJa+8qB6AISEZkkim3d9N3wclN2eqqZ3RZ5VBEbTKZpHmeS2LH/CIeGhiOMSkSk8hVbJ3G2ux/ITrj7foIWThVtMJGmdYyXmyBoBguwYafqJUSkthSbJGLhzXQAmNk0xnHJaqIMJFO0jLF1E8CyOVMAeHbHwahCEhGZFIo9wX8ReNzMfhxOvx/4XLQhRW8wOb6SRHd7I3M6m1m740B0QYmITALFtm76FzNbDfxWOOt97r4++rCiNZAYX0kCYPn8TtZuPxBNQCIik0TR/16HSaHiE0NWOuMkUplxlSQAls/t5KfP7qL3cILu9saIohMRqWxVP57EYDLot6llHK2bIChJADyjprAiUkNqIEmkAcZ8M13WslOmUBcznlG9hIjUkHElCTObbWYVfe1lYBw9wOZqbqhjycx23VQnIjVlvCWJ7wEbzewLUQRTClGVJOBo5XVa3YaLSI0YV5Jw98uARcC3owkneiMliXHWSQC8aeE0DidSrNdNdSJSI8ZdJ+GB56MIphQGh4OSxHi65chauagLgCde6hv3tkREJoOirsGYWRPwx8CbAQceBW7JjjExHma2BPhRzqxFwF+7+5fHs93BRJAkxtpVeK4ZHU0s6m7l8Zf6+OhbF417eyIixXIPmvX3J1L0D6XoT6QYSATPwes0/Ylh+hPpYP5Qiv5k6ujr8OpKoYo9c/4LcBj4Wjj9QYJ6ifcXuZ3juPsmYDmAmdUBrwJ3jXe7AxE1gc1asaiL/1i7k1Q6M6aR7kSkNrg7Q8MZBpMpBpNpjgyng+dkmiPDwbzsdPAczhsO10mmGQhP7sGJ/2hCSBVYL9raUEdrY5y2xjhtTXFaG+LMndpS1OcoNkksc/elOdMPmVkpbqy7FNji7uMe0GhwpE4imi6mVizq4o4nt7N+1yHOntsZyTZFpLTcneG0M5RKkxjOkEilGQqfE6kMQ8PBc3ZZYrRl2dfHvH9oOBOe2FMcGc456Q+n8SLbuDTX19HSUEdzQ/Y5TltjHd3tjbQ2xmlvjNMaPtrDk/7I68Zg3bbGelob62htiBOL2aj7+daHCo+p2DPnU2a2wt2fADCzNwGri9xGIa4GfjDaAjO7AbgBYP78+Sfd0EC2ddM4m8BmrVg4DYDHt/QpSciEcHcyDhl3Mu74yOvg2TNHl2W8gPVzl4fvhTzrZMawzZH1j27/uPWznyvcfjrjJNMZhtMZUmlnOJ0ZmR5OOcOZDMNpZziVGVmWXS+Y9vC9wXpH35thOBMsK/aEncsMmuJ1NNbHaIzHaKqvozEeozEePHc0xZnV0UhLQzw4wYcn+6aR1/GcE38w3dJQR3P90YTQFK/Le1Ivp2KTxHnAY2a2PZyeB7xgZs8R1GGfPd6AzKwBeBdw02jL3f1W4FaAnp6ek37tg8kU8ZjRENGloRkdTZw5q52HNu3hYxctjmSbcnKpdIaBZDoseqcYSKYZTKQYSqVJhv/lJVPBySOZynmE04nw5JJxJ51x0uHJK51x0uHJKp3xnHmMzEuPnMxGO2HmnBD92BNicJIsav1wee5JtdbUxYz6OqO+LhY+gtcN4XQ8Z7qxPkZbU5x4LEZDPPc9seO2cfTEHqMx5yTfVB+e7Otjr0sEucvq6wyzyjuBT4SCkoSZnQbMBK48ZtE8YDcQ5Wg8vw085e6vRbGx7IBDUX7Bl541g68//BIHB4eZ0lIf2XZrRSqdobc/wa6DQ7x2cIi9/Qn2DQyzfzDJgcEk+wezr4dHrsEmUpkx7csMGupiNMSDk0osZtSZURczYjGoM3v9vJFlRp0xMq8hHiMWrhszgtcGZrnThuUsC6ZzlsdOvj4nev/I/nKX56wfK2b97PJC1jm6/cK2OcpnzL6X0deJxRg5oddV4H/TtazQksSXgZuOrSMwsw7gS+7+zghjuoY8l5rGYrwDDo3m0rNm8o8PbeEXL+zh3cvnRLrtauDuvHYowda9A2zrG2Db3gFe7htk18Ej7D40RO/hxKj/Ibc3xulsrWdqSwOdLQ0s6GqlvSmodGtpiAfXWRuPXnttaYjTVF83kgQa40eTQUP4Oh6r3f8ARaJQ6Nlzprs/d+xMd3/OzBZEFYyZtQKXAx+LapvjHXBoNMvndjK9rYFVG5QkBpMpNuw6zIZdh1i/6xAbdh1i0+7DI3e6Q/Cf/LxpzZzS2cySWe3M6mhi1pRmZk1pZFZHM93tjXS21FOv1mIiFafQJNF5gmXNEcQBgLsPAF1RbQ/GP+DQaGIx45IlM7j/+d0MpzM1dXLbdfAIq7ftZ/W2fax+eT8bdh0aKRW0N8VZOruDD/TMY3F3Kwumt7Kgq5VTOpt1CUFkkir07LnazD7q7t/InWlmHwHWRB9WdAYSqcjukch12dKZ/HjNDh7b0sdFZ3RHvv1KMZhM8djmPn7xwh4efqGXV/YdAYKmeufO7+QTl5zGb8zt5KzZ7czpbNalHZEqU2iS+BRwl5n9LkeTQg/QALy3BHFFZjCZZnpbQ+TbvXhJN1Oa67nrqR1VlyT2HB7iP5/dxaqNe3jypX0k0xma6+v4zdO6uO7ChfQsmMpZsztqqgQlUqsKShJhS6MLzewSYFk4+6fu/l8liywig8kULY3F3WFYiMZ4HW8/ezZ3PfUq/YkUbRF0+1FOBweH+dm6XdzzzE6eeKmPjMPi7lb+YOWpXLxkBucvnEpjPPoSmYhUtmLHuH4IeKhEsZREUCdRmpPb+86dwx1Pbue+dbu56ry5JdlHKWUyzq+27OX2J7azauNrDKedBV0tfOKS03jnOadw+sz2cocoImU2uf/9LUBQJ1Gaj3neqVOZP62FO9e8MqmSxL6BJHeueYU7ntzOtr5BprU2cO3KBbx7+RyWzelQvYKIjKjqJOHuQUki4iawWWbGNRfM5/P3bWT9zkMsPaWjJPuJyuY9/fzzw1u4+5mdJFMZLlgwjU9ffgZXLpulS0kiMqqqThLJdIZUxktWkgD44AXz+eqqF/nWo1v54gfOKdl+xuPp7fv5+sNbeGD9azTGY/xOzzx+b8WpLJmly0kicmJVnSSOjAxdWrr/kqe01POBnrnc8evt/PmVS5jZ0VSyfRXrsS17+dqqzTz+Uh8dTXE+eclpXHvhArraKnpYchGpIFWdJLI9wEZ9M92xPvzmhXzviZe55RdbuPldbyjpvgrx1Pb9fPGBTfxqcx8zOxr5q7efxdUXzJ/0LbBEZOJV9VkjO5ZE1N1yHOvUrlZ+5/x5fP+Jl7n2wgUsnN5a0v3ls37nIb74wCZWbdxDV2sDf/X2s/i9FafSVK/6BhEZm6pOEv0RDzh0Ip++/AzuXruTz/9sI1///fNKvr9cW3r7+dKDL3Dvs7tob4rzP684g+t+c2EkQ7aKSG2r6rPIQITjW5/MjPYm/uiixXzxwRd44PndXPGGWSXf56sHjvCVn7/AnWt20FRfxycuOY2PvmWRui8XkchUdZIYKUmU+HJT1scuWsx9z+/mMz95juXzO5nRXppK7L39Cf7xoc3c/kQw9tOHLlzIH1+ymOmqkBaRiFV1khgIk0R748T8Z90Qj/GVq5fz9q8+yifveJrvfviCSOsDDg0N841HXuJbj25laDjN+8+bx59cdjpzOiPriFdE5HWqOklMdEkC4LQZ7fyfq87mxh+u5VM/XMv/++C5xMfZEd7BI8N87/FtfOOXWzl4ZJi3nz2bP738DBZ3t0UUtYjI6GokSUzsx3z38jn09Sf523vXc/13V/PVq88dUz3B7oNDfPtXW7n9ye30J1JcsqSbP7tiCcvmTClB1CIix6vqJDGQSBGPGY3xie/S+sNvXkhzQx1/ffc63v61X/K/37GUK5bOPGm/SIlUmlUb9vDj1a/w8Au9ALzznFP42FsXV3y3HyJSfao+SbQ1xcvWYd01F8znjJltfObfnuNj31vDmbPaec+5c+g5dSrzu1porq+jP5Fi54Eh1u86xBMv9fHoi3s5eGSYWR1N/OFFi7n6/PnM74q+q3MRkUJUdZI4nEhNyD0SJ3LeqdP42Y1v4cdrdvDD/36Fv//Zxrzrzmhv5PKlM3nH2bN5y+ndGvJTRMquqpPEQIUMBhSvi3HNBfO55oL57D44xPM7D7LzwBGGhjO0NsaZ0d7IWad0cMqUJnXTLSIVpfxn0BIaSKRpa6qsjzhrShOzplROJ4AiIidS1YMUH06k1DWFiMg4VHWSCC43qXM7EZGxqvokUe6KaxGRyayqk0T/UKri6iRERCaTqk0S7s5AsjJaN4mITFZVmySODKfJ+MR3ySEiUk0qKkmYWaeZ3WlmG81sg5mtHOu2sv02qSQhIjJ2lXYG/Qpwn7tfZWYNwJj7o+gfUpIQERmvijmDmtkU4K3AhwDcPQkkx7q9iRyVTkSkWlXS5aaFQC/wbTN72sy+aWatY91YOcaSEBGpNpWUJOLAG4Fb3P1cYAD4zLErmdkNZrbazFb39vbm3Vj/BI9KJyJSjSopSewAdrj7k+H0nQRJ43Xc/VZ373H3nu7u7rwbG1BJQkRk3ComSbj7buAVM1sSzroUWD/W7al1k4jI+FXaGfSTwO1hy6aXgOvGuqGBMg1dKiJSTSrqDOrua4GeKLbVn0hhBi0NutwkIjJWFXO5KWr9iRRtDeUbulREpBpUbZIY0FgSIiLjVrVJoj+hHmBFRMaripNEWiUJEZFxqtokoVHpRETGr6qThEalExEZn6pNEvsGkkxtaSh3GCIik1pVJolMxtk3kKSrTUlCRGQ8qjJJHBoaJpVxutoayx2KiMikVpVJom8gGIZiukoSIiLjUp1Joj9IEl2tKkmIiIxHlSaJBADTWlWSEBEZj6pMEnt1uUlEJBJVmSSyJYmpKkmIiIxLVSaJfQNJOlvqqa+ryo8nIjJhqvIs2tefVH2EiEgEqjJJ7O1PMF0tm0RExq0qk0Sf7rYWEYlEdSaJ/oSShIhIBKouSaTSGfYPDutGOhGRCFRdktg/OAygkoSISASqLkn0DQT3SKgkISIyftWXJLL9NqkkISIyblWXJPaGd1urSw4RkfGruiSRLUlM0+UmEZFxq7oksW7nQaa1NtDZXF/uUEREJr2qShLuzhNb+li5qItYzModjojIpFdVSeLlvkF2Hhxi5eKucociIlIV4uUOIJeZbQMOA2kg5e49xbz/sS19AFyoJCEiEomKShKhS9x971je+NiWvczqaGLh9NaoYxIRqUlVc7npwGCSx7f0sXJxF2aqjxARiYK5e7ljGGFmW4H9gAP/7O63jrLODcAN4eQyYN3ERViRpgNjKnlVER2DgI6DjgEUdgxOdffuQjZWaUlijru/amYzgAeBT7r7IydYf3Wx9RbVRsdAxyBLx0HHAKI/BhV1ucndXw2f9wB3AReUNyIRkdpWMUnCzFrNrD37GrgCXUoSESmrSmrdNBO4K6x0jgN3uPt9J3nPcXUWNUjHQMcgS8dBxwAiPgYVVSchIiKVpWIuN4mISOVRkhARkbwmZZIwsyvNbJOZbTazz5Q7nlIxs3lm9pCZrTez583sxnD+NDN70MxeDJ+nhvPNzL4aHpdnzeyN5f0E0TGzOjN72szuDacXmtmT4Wf9kZk1hPMbw+nN4fIFZQ08QmbWaWZ3mtlGM9tgZitr7bdgZp8O/xbWmdkPzKypFn4LZnabme0xs3U584r+7s3s2nD9F83s2kL2PemShJnVAf8I/DawFLjGzJaWN6qSSQF/5u5LgRXAx8PP+hlglbufDqwKpyE4JqeHjxuAWyY+5JK5EdiQM/154EvufhrBDZjXh/OvB/aH878UrlctvgLc5+5nAucQHI+a+S2Y2RzgT4Aed18G1AFXUxu/he8AVx4zr6jv3symAZ8F3kRwe8Fns4nlhNx9Uj2AlcD9OdM3ATeVO64J+ux3A5cDm4DZ4bzZwKbw9T8D1+SsP7LeZH4Ac8M/gt8C7gWM4I7S+LG/CeB+YGX4Oh6uZ+X+DBEcgynA1mM/Sy39FoA5wCvAtPC7vRf4H7XyWwAWAOvG+t0D1xD0ZMFo6+V7TLqSBEd/KFk7wnlVLSwqnws8Ccx0913hot0EzYeheo/Nl4E/BzLhdBdwwN1T4XTu5xw5BuHyg+H6k91CoBf4dnjZ7Zvh/UQ181vw4GbbLwDbgV0E3+0aau+3kFXsdz+m38RkTBI1x8zagH8DPuXuh3KXefAvQdW2YzazdwB73H1NuWMpszjwRuAWdz8XGODo5QWgJn4LU4F3EyTMU4BWjr8EU5NK+d1PxiTxKjAvZ3puOK8qmVk9QYK43d1/Es5+zcxmh8tnA3vC+dV4bH4TeFc41sgPCS45fQXoNLPszaC5n3PkGITLpwB9ExlwiewAdrj7k+H0nQRJo5Z+C5cBW929192HgZ8Q/D5q7beQVex3P6bfxGRMEv8NnB62aGggqLi6p8wxlYQFt59/C9jg7v+Qs+geINsy4VqCuors/D8IWzesAA7mFEcnJXe/yd3nuvsCgu/6v9z9d4GHgKvC1Y49Btljc1W4/qT/79rddwOvmNmScNalwHpq6LdAcJlphZm1hH8b2WNQU7+FHMV+9/cDV5jZ1LBUdkU478TKXRkzxgqctwEvAFuAvyx3PCX8nG8mKEI+C6wNH28juK66CngR+DkwLVzfCFp+bQGeI2gFUvbPEeHxuBi4N3y9CPg1sBn4MdAYzm8KpzeHyxeVO+4IP/9yYHX4e/h3YGqt/RaAvwE2EvTr9j2gsRZ+C8APCOphhglKldeP5bsHPhwej83AdYXsW91yiIhIXpPxcpOIiEwQJQkREclLSUJERPJSkhARkbyUJEREJC8lCZnUzKzLzNaGj91m9mrOdEO548tlZheb2YUTtK9tZjY9fH2emW01s3MnYt9SXSpp+FKRorl7H8H9A5jZzUC/u3+hXPGYWdyP9iN0rIuBfuCxiLZXyPvPJrg7+3fc/emxbkdql0oSUnXC/5wfNrM1ZnZ/TtcFvzCzL5nZagvGYzjfzH4S9q3/d+E6CywYr+H2cJ07zaylgO1+2cxWAzea2TvD8QueNrOfm9nMsIPGPwQ+HZZy3mJm3zGzq3Li7g+fLzazX5rZPcB6C8bS+L9m9t/h+AAfK/BQnEVw093vu/uvIzm4UnOUJKTaGPA14Cp3Pw+4DfhczvKku/cAXyfoxuDjwDLgQ2aW7SF0CfBP7n4WcAj447APrRNtt8Hde9z9i8CjwAoPOuL7IfDn7r4t3OeX3H25u//yJJ/jjcCN7n4Gwd21B939fOB84KNmtrCAY3E38Al3f7SAdUVGpctNUm0aCU76Dwbd+1BH0J1BVrafr+eA5z3sz8jMXiLo/OwA8Iq7/ypc7/sEA93cd5Lt/ijn9VzgR2FJo4FgHIhi/drds++7Ajg7p9QxhWBAmZNt9+fAR8zsfndPjyEGESUJqTpGcPJfmWd5InzO5LzOTmf/Ho7tq8YL2O5AzuuvAf/g7veY2cXAzXnekyIszZtZjCChjLY9Az7p7ifvjO31PkFQevknoNBLVCKvo8tNUm0SQLeZrYSgq3Uze0OR25iffT/wQYLLR5uK2O4UjnbBnDuO8GGgPWd6G3Be+PpdQH2e7d0P/FF4yQszO8OCAYcws40n+ByZMP4zzexvT7CeSF5KElJtMgTdQn/ezJ4h6Dm32GanmwjGE99A0NPqLe6eLGK7NwM/NrM1BENmZv0H8N5sxTXwDeCicHsreX3pIdc3CbrEfsrM1hEMOxkPm7jaiT6Iuw8RJKB3mdnHT/yxRY6nXmBFcoStkO5192XljuVkLBi1b5G7f7XcsUj1Up2EyCTl7veWOwapfipJiIhIXqqTEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8/j95GDoyqNqccwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q = DiatomicPartitionFunction(mol, energy_surface, vibrational_structure)\n",
    "\n",
    "P = 101350  # Pa\n",
    "temps = np.arange(10, 1050, 5)  # K\n",
    "\n",
    "mol.spins = [1 / 2, 1 / 2]\n",
    "\n",
    "td = Thermodynamics(Q, pressure=101350)\n",
    "td.set_pressure(101350)\n",
    "temps = np.arange(10, 1500, 5)\n",
    "ymin = 5\n",
    "ymax = 11\n",
    "\n",
    "plt.plot(temps, td.constant_pressure_heat_capacity(temps) / const.CAL_TO_J)\n",
    "plt.xlim(0, 1025)\n",
    "plt.ylim(ymin, ymax)\n",
    "plt.xlabel(\"Temperature, K\")\n",
    "plt.ylabel(\"Cp, cal mol$^{-1}$ K$^{-1}$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we demonstrate how to access particular components (the rotational part) of the partition function, which in the H2 case we can further split to para-hydrogen and ortho-hydrogen   components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq = Q.get_partition(part=\"rot\", split=\"eq\")\n",
    "para = Q.get_partition(part=\"rot\", split=\"para\")\n",
    "ortho = Q.get_partition(part=\"rot\", split=\"ortho\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now plot the constant volume heat capacity (of the rotational part) demonstrating how we can call directly the functions in the 'thermodynamics' module, providing a callable object for the partition function (or in this case its rotational component). Note that in the plot we normalize the plot dividing by the universal gas constant R (Avogadro's number times Boltzmann's constant) and we use crossed to compare with experimental data found in literature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# REFERENCE DATA from literature\n",
    "df_brink_T = [80.913535, 135.240157, 176.633783, 219.808499, 246.226899]\n",
    "df_brink_Cv = [0.118605, 0.469925, 0.711510, 0.833597, 0.895701]\n",
    "\n",
    "df_eucken_T = [\n",
    "    25.120525,\n",
    "    30.162485,\n",
    "    36.048121,\n",
    "    41.920364,\n",
    "    56.195875,\n",
    "    62.484934,\n",
    "    72.148692,\n",
    "    73.805910,\n",
    "    73.804236,\n",
    "    92.214423,\n",
    "    180.031917,\n",
    "    230.300866,\n",
    "]\n",
    "df_eucken_Cv = [\n",
    "    0.012287,\n",
    "    0.012354,\n",
    "    0.008448,\n",
    "    0.020478,\n",
    "    0.032620,\n",
    "    0.048640,\n",
    "    0.048768,\n",
    "    0.076678,\n",
    "    0.078670,\n",
    "    0.170548,\n",
    "    0.667731,\n",
    "    0.847681,\n",
    "]\n",
    "\n",
    "df_gia_T = [\n",
    "    190.919338,\n",
    "    195.951254,\n",
    "    202.652107,\n",
    "    204.292585,\n",
    "    209.322828,\n",
    "    225.300754,\n",
    "    234.514217,\n",
    "    243.747768,\n",
    "]\n",
    "df_gia_Cv = [0.711700, 0.723719, 0.749704, 0.797535, 0.811546, 0.797814, 0.833793, 0.845868]\n",
    "\n",
    "df_parting_T = [80.101665, 86.358919, 185.914204, 239.927797]\n",
    "df_parting_Cv = [0.084730, 0.138598, 0.667809, 0.891634]\n",
    "\n",
    "df_ce_T = [\n",
    "    80.669344,\n",
    "    135.550569,\n",
    "    145.464190,\n",
    "    165.301153,\n",
    "    182.144856,\n",
    "    203.372528,\n",
    "    237.993108,\n",
    "    268.696642,\n",
    "    294.095771,\n",
    "    308.872014,\n",
    "]\n",
    "df_ce_Cv = [\n",
    "    0.103048,\n",
    "    0.467344,\n",
    "    0.541364,\n",
    "    0.647315,\n",
    "    0.714078,\n",
    "    0.798258,\n",
    "    0.891147,\n",
    "    0.944848,\n",
    "    0.966618,\n",
    "    0.985486,\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAABfcklEQVR4nO3deXxNV9fA8d/OHIl5iprHIiRBzLOi1PRQiqLSQYtSOlGtt0KVFi0tWqVqqpZWTTVWK4IaY4qpIpEYgwgRIfPd7x/7JoLMcnNvkv39fM6TO5ycs3LrycreZ5+1hJQSTdM0TbM0VuYOQNM0TdNSoxOUpmmaZpF0gtI0TdMskk5QmqZpmkXSCUrTNE2zSDbmDiCrSpUqJatUqWLuMDRN07QMHDly5JaUsnR2vz/PJagqVarg5+dn7jA0TdO0DAghLj7N9+spPk3TNM0i6QSlaZqmWSSdoDRN0zSLpBOUpmmaZpF0gtI0TdMskk5QmqZpmkXSCUrTNE2zSDpBaZqmaRZJJyhN0zTNIukEpWmaxbl+/ToDBgygevXqNGrUiBdeeIGAgACTnGvXrl3s27cv1feWLl1K6dKl8fDwSN7OnDmTrfNs3LiRL774AgBvb29mzZoFwKeffsrff/8NqEo5t27dytKx8rM8V+pI07T8TUpJ7969GTp0KKtWrQLgxIkT3Lhxg1q1amXrmAkJCdjYpP7rbteuXTg7O9OiRYtU3+/fvz/z5s3L1nlT6tmzJz179nzi9SlTpmTpOAkJCWkeK7/RCUrTtDSNHTuW48eP5+gxPTw8mDNnTprv+/j4YGtry/Dhw5Nfc3d3B2DAgAEMGTKEbt26AeDl5UX37t3p27fvE8dZunQpa9euJSoqisTERNatW8drr73GhQsXKFSoEAsXLqRIkSIsWLAAa2trfv75Z+bOnUvr1q0z/BmklIwePZodO3ZQsWJF7OzseO211+jbt29yvdBSpUrh5+fHBx98wK5du1i6dCl+fn5PJLvHf4YZM2awdetWHB0d+eWXX6hRowZeXl44ODhw7NgxWrZsiZubW/KxHv9+Z2dnoqKi2LVrF5MmTaJYsWKcPHmSl156ifr16/PNN98QHR3N+vXrqV69eoY/qznpBKVpmkU5deoUjRo1SvW9/v3789tvv9GtWzfi4uL4559/+P7779M81tGjR/H396dEiRKMHj2aBg0asH79enbu3Mkrr7zC8ePHGT58OM7OznzwwQepHmP16tXs3bs3+fn+/fvZunUr586d48yZM9y4cYO6devy2muvPd0PblS0aFFOnjzJ8uXLGTt2LJs2bQLgypUr7Nu3D2tra5YuXZqpY504cYKzZ89SokQJqlWrxhtvvMGhQ4f45ptvmDt3brp/KFgCnaA0TUuTpf0C69q1K2PGjCE2NpZt27bRpk0bHB0d09y/U6dOlChRAoC9e/fyxx9/ANChQwfCw8OJjIzM8JypTfHt3r2bgQMHYm1tzTPPPEOHDh2e4qd61MCBA5O/vvvuu8mv9+vXD2tr6ywdq3HjxpQrVw6A6tWr07lzZwDq16+Pj49PDkVsOnqRhKZpFsXV1ZUjR46k+p6DgwPt2rVj+/btrF69mv79+6d7LCcnJ1OEmC4bGxsMBgMAMTExWf5+IUSqj9P6WVKez2AwEBcXl/yevb198mMrK6vk51ZWViQkJGQ5ttymE5SmaRalQ4cOxMbGsnDhwuTX/P392bNnD6BGNEuWLGHPnj106dIl08dt3bo1K1euBNTCiFKlSlGkSBEKFy7MvXv3shRjmzZtWL16NYmJiYSGhj4yGqlSpUpygk0asWXF6tWrk782b948w/1Tnm/jxo3Ex8dn+ZyWSicoTdMsihCCdevW8ffff1O9enVcXV2ZMGECLi4uAHTu3BlfX186duyInZ1dpo/r7e3NkSNHcHNz46OPPmLZsmUA9OjRg3Xr1uHh4ZGcBFNavXr1I8vM9+3bR+/evalZsyZ169bllVdeeSSRTJo0iTFjxuDp6ZnlKTmAO3fu4ObmxjfffMPs2bMz3H/YsGH4+vri7u7O/v37zTJqNBUhpTR3DFni6ekpdUddTdMsSXqrCQsyIcQRKaVndr9fj6A0TdM0i6RX8Wmalqdt376d8ePHP/Ja1apVWbduXa7FkNll31rWmGyKTwhREVgOlAUksFBK+c1j+7QDNgDBxpfWSinTva1aT/FpmqblDU87xWfKEVQC8L6U8qgQojBwRAixQ0r5eCGrPVLK7iaMQ9M0TcuDTHYNSkoZKqU8anx8DzgLlDfV+TRN07T8JVcWSQghqgANgIOpvN1cCHFCCLFVCOGaxve/KYTwE0L4hYWFmTJUTdM0zUKYPEEJIZyBP4CxUsrH64ocBSpLKd2BucD61I4hpVwopfSUUnqWLl3apPFqmqZplsGkCUoIYYtKTiullGsff19KGSmljDI+3gLYCiFKmTImTdMsnyX2g6pbty6LFi0ySQxa6kyWoIQqIrUYOCul/DqNfVyM+yGEaGKMJ9xUMWmaZvmS+kG1a9eOoKAgjhw5wvTp07lx40a2j5le3bn0EhSo0krHjx9n165dfPzxx5mOIzExMctxao8y5Sq+lsAQ4KQQ4rjxtY+BSgBSygVAX2CEECIBiAYGyLxW2kLT8rGxYyGH20Hh4QHpFUm31H5QZcqUoXr16ly8eBFvb28OHz5MdHQ0ffv2ZfLkyYCqi9e/f3927NjBuHHjuHfvHgsXLiQuLo4aNWqwYsUKChUqlL0PrgAyWYKSUu4FRAb7zAOevlWlpmn5hqX1g0py4cIFLly4QI0aNfj8888pUaIEiYmJPPfcc/j7++Pm5gZAyZIlOXr0KADh4eEMGzYMgIkTJ7J48WJGjx6dnY+lQNKVJDRNS5OFtYMySz+opIaF9vb2/PDDD5QoUYIFCxawcOFCEhISCA0N5cyZM8kJKmULkFOnTjFx4kQiIiKIiori+eeff5ofv8DRCUrTNIvi6urKmjVrUn3v8X5QAwYMSPdYOVHZ+/GGhcHBwcyaNYvDhw9TvHhxvLy8Hun7lPKcXl5erF+/Hnd3d5YuXcquXbueOp6CRBeL1TTNolh6P6jIyEicnJwoWrQoN27cYOvWrWnue+/ePcqVK0d8fHzyubXM0wlK0zSLYmn9oB7n7u5OgwYNqF27Ni+//DItW7ZMc9/PPvuMpk2b0rJlS2rXrp3pWDVF94PSNE3TTEL3g9I0TdPyJb1IQtO0PM0S+kFppqGn+DRN0zST0FN8mqZpWr6kE5SmaZpmkXSC0jRN0yySTlCapmmaRdIJStM0i2Mp/aAA1q9fj5ubG3Xq1KF+/fqsX78+3X3PnDmT/Lxdu3boRV3ZpxOUpmkWxZL6QZ04cYIPPviADRs2cPbsWTZu3MgHH3yAv79/qud4PEFpT0cvM9c0LW1maAi1c+dOvL292b179xPvmaIfVLNmzbC2tqZ06dJP9IMaMmQI7du357XXXkt+bfHixezatYsVK1bQrl07PDw82Lt3L7179+arr76iaNGiFC1alD/++IPXX3+dpk2b4uPjQ0REBIsXL6Z169bExMQwYsQI/Pz8sLGx4euvv6Z9+/bZ/0wt1NMuM9c36mqaZlEsqR/U6dOnn3jd09OT+fPnJz+Pi4tLnsY7f/78EwkzISGBQ4cOsWXLFiZPnszff//N/PnzEUJw8uRJ/vvvPzp37kxAQAAODg5Z+qzyO52gNE1Lm4U1hDJHP6iMpOz/lJo+ffoA0KhRI0JCQpJjSWpcWLt2bSpXrkxAQEByTylN0degNE2zKK6urhw5ciTV9x7vB5VRcnjaflB169Z9IpYjR47g6uqa6XPY29sDYG1tne61MO1JOkFpmmZRLKkf1AcffMD06dOTRz4hISFMmzaN999/P9X9M9tbKmUsAQEBXLp0iWeffTbTP0tBoROUpmkWxZL6QXl4ePDll1/So0cPateuTY8ePZgxYwYeHh6pnmPAgAHMnDmTBg0aEBQUlGYsI0eOxGAwUL9+ffr378/SpUuTR1raQ3oVn6ZpmmYSulispmmali/pVXyapuVpuh9U/qWn+DRN0zST0FN8mqZpWr6kE5SmaZpmkXSC0jRN0yySTlCapmmaRdIJ6imFh4fz2muvsWrVKl3GRNNyiKX0g/rvv/9o3rw59vb2zJo1K81jvP7667i7u+Pm5kbfvn2JiooySaw56fGf28vLizVr1pgxoifpBPWUPvnkE5YsWcLAgQOpXr06c+bMyVSpE03TUmdJ/aBKlCjBt99+m2ql85Rmz57NiRMn8Pf3p1KlSsybNy/bsaZGSonBYMix4yUkJGTYqNES6PugnoK/vz+LFi1i1KhRdO7cmVmzZvHuu+/i7e3N+PHjmTBhgrlD1LSnMnbbWI5fP56jx/Rw8WBOlzlpvu/j44OtrS3Dhw9Pfs3d3R0wTT+oBQsWYG1tzc8///xEP6gyZcpQpkwZNm/enO7PVKRIEUAlkujoaIQQT+zj7e1NUFAQgYGB3Lp1i3HjxjFs2DCioqLo1asXd+7cIT4+nqlTp9KrVy9CQkJ4/vnnadq0KUeOHGHLli188cUXHD58mOjoaPr27cvkyZOfOE9ISAivvfYat27donTp0ixZsoRKlSrh5eWFg4MDx44do3z58uzbt++Rnxtg9+7dfP3111y/fp0ZM2bQt29fpJSMGzeOrVu3IoRg4sSJGRbpzSkmG0EJISoKIXyEEGeEEKeFEGNS2UcIIb4VQgQKIfyFEA1NFU9Ok1IyduxYihcvzpQpU+jRowe+vr4cPHiQli1b8vHHH+Pr62vuMDUtz8lMPygguR9UUrJKzdGjR1mzZg2+vr5MmjSJBg0a4O/vz7Rp03jllVeoUqUKw4cP59133+X48eOPJKeMvPDCC1y7di35+auvvoqLiwv//fdfciuNx/n7+7Nz507279/PlClTuHbtGg4ODqxbt46jR4/i4+PD+++/T9L9qefPn2fkyJGcPn2aypUr8/nnn+Pn54e/vz++vr7JnX0//fRTNm7cCMDo0aMZOnQo/v7+DBo0iHfeeSf5/FeuXGHfvn2sXbs21Z87NDSUvXv3smnTJj766CMA1q5dy/Hjxzlx4gR///03H374IaGhoZn+nJ6GKUdQCcD7UsqjQojCwBEhxA4pZcp+yF2BmsatKfC98avFW79+PT4+PsyfP5/ixYsnv96kSRPWrFlDzZo1GTduHAcOHEj1rylNywvSG+mYgyX1g9qyZcsjz5csWUJiYiKjR49m9erVvPrqq098T69evXB0dMTR0ZH27dtz6NAhunXrxscff8zu3buxsrLi6tWrydOZlStXplmzZsnf/9tvv7Fw4UISEhIIDQ3lzJkzuLm5MWXKlOR99u/fz9q1awHVEXjcuHHJ7/Xr1w9ra+s0f6b//e9/WFlZUbdu3eQY9u7dy8CBA7G2tqZs2bK0bduWw4cP07Nnz2x8alljshGUlDJUSnnU+PgecBYo/9huvYDlUjkAFBNClDNVTDklJiaG999/n3r16vHmm28+8b6joyOfffYZhw4dsriLjppm6SypH1RWWVtbM2DAgORE+LjH/1gVQrBy5UrCwsI4cuQIx48fp2zZssTExACPxh8cHMysWbP4559/8Pf3p1u3bsn7ZVZme1cBWEKVoVxZJCGEqAI0AA4+9lZ54HKK51d4MokhhHhTCOEnhPALCwszWZyZNWfOHIKDg5k9ezY2NqkPQl955RXq1avHxx9/THx8fC5HqGl5lyX1g8oMKSWBgYHJjzdu3Ejt2rVT3XfDhg3ExMQQHh7Orl27aNy4MXfv3qVMmTLY2tri4+PDxYsXU/3eyMhInJycKFq0KDdu3GDr1q2p7teiRQtWrVoFwMqVK9OctsxK76rVq1eTmJhIWFgYu3fvpkmTJhl+X46QUpp0A5yBI0CfVN7bBLRK8fwfwDO94zVq1Eia07Vr16Szs7Ps2bNnhvtu2rRJAnLevHm5EJmm5R9Xr16V/fr1k9WqVZN169aVL7zwggwICJBSShkXFyeLFy8uvby80j3GkiVL5Ntvv538PDw8XPbq1UvWr19fNm3aVJ44cUJKKeW5c+dk/fr1pbu7u9y9e/cjxwgNDZXly5eXhQsXlkWLFpXly5eXd+/elVJK2bVrV3n16lWZmJgoW7RoIevVqyddXV3lyy+/nLxPSpMmTZJDhgyRzZo1kzVq1JALFy6UUkoZFhYmmzVrJuvVqye9vLxk7dq1ZXBwsAwODpaurq6PHGPo0KGyZs2askOHDrJ3795yyZIlUkop/+///k9u2LBBSillSEiIbN++vaxfv77s0KGDvHjxYvL3/v7778nHevznfvx9JycnKaWUBoNBfvDBB9LV1VXWq1dPrlq1Kt3PPSXATz5F/jBpsVghhK0xCW2XUn6dyvs/ALuklL8an58D2kkp07wCZ+5isa+//jorVqzg9OnT1KxZM919pZS0b9+eM2fOEBQUROHChXMpSk3TLI23tzfOzs4ZLlnPTyy2WKxQk62LgbOpJSejjcArxtV8zYC76SUnc4uPj+fXX3/l1VdfzTA5gZpfnjFjBmFhYene5KdpmqY9yZSr+FoCQ4CTQojjxtc+BioBSCkXAFuAF4BA4AHw5LIXC3Ly5Emio6Np3759pr+nSZMm9OvXj6+++ooRI0Ykt63WNC1n5JV+UN7e3uYOIc/R/aCyYP78+YwaNYqQkBAqV66c6e8LDAykTp06jBw5km+++caEEWqaplkOi53iy48OHDiAi4sLlSpVytL31ahRg5deeolly5bx4MEDE0WnaZqWv+gElQX79++nefPm2brxdtiwYdy9e1ffF6VpmpZJOkFlUlhYGEFBQY/c1Z0Vbdu2pWbNmixatCiHI9M0TcufdILKpAMHDgDQvHnzbH2/EII33niDvXv3cvbs2ZwMTdO0p/Tpp5/y999/m/QcS5cufaR2X1oy0/YiJCSEevXqZbjPL7/8kqUYLY1OUJl04MABbGxs0iximRlDhw7FxsaGH3/8MQcj0zTL4L3L29whZEtiYiJTpkyhY8eOJj1PZhNUTtEJqgDZv38/7u7uFCpUKNvHKFu2LL169WLZsmXExsbmYHSaZn6TfZ9s/ZBdP//8M02aNMHDw4O33nqLxMREDh8+jJubGzExMdy/fx9XV1dOnTrFrl27aNOmDd26dePZZ59l+PDhyb2T/vrrL5o3b07Dhg3p169fciPBKlWqMH78eBo2bMjvv//+yKilSpUqTJgwAQ8PDzw9PTl69CjPP/881atXZ8GCBckxzpw5k8aNG+Pm5sakSZMAlRTq1KnDsGHDcHV1pXPnzkRHR7NmzRr8/PwYNGgQHh4eREdHM2XKFBo3bpxc0zOjFdVHjhzB3d0dd3d35s+fn/x6SEgIrVu3pmHDhjRs2DC5x9NHH33Enj178PDwYPbs2WnuZ9GepgyFOTZzlDpKSEiQTk5Oj5RNya5t27ZJIEvlQjQtL8CbHDnOmTNnZPfu3WVcXJyUUsoRI0bIZcuWSSml/OSTT+T7778vR44cKadNmyallNLHx0fa29vLoKAgmZCQIDt27Ch///13GRYWJlu3bi2joqKklFJ+8cUXcvLkyVJKKStXriy//PLL5HOmLPNTuXJl+d1330kppRw7dqysX7++jIyMlDdv3pRlypSRUkq5fft2OWzYMGkwGGRiYqLs1q2b9PX1lcHBwdLa2loeO3ZMSillv3795IoVK6SUUrZt21YePnw4+Zzh4eHJjwcPHiw3btz4RCwp1a9fX/r6+kopZXLpISmlvH//voyOjpZSShkQECCTfkf6+PjIbt26JX9/WvuZEk9Z6kg3LMyEU6dOcf/+/Wxff0qpU6dOVK5cmUWLFuVa0y9NMxXvXd6PjJzEZLXCdVLbSXi3887WMf/55x+OHDlC48aNAYiOjqZMmTKAulbUuHFjHBwc+Pbbb5O/p0mTJlSrVg2AgQMHsnfvXhwcHDhz5gwtW7YEVP+olP8fTu//f0mtJOrXr09UVBSFCxemcOHC2NvbExERwV9//cVff/1FgwYNAIiKiuL8+fNUqlSJqlWr4uHhAUCjRo0ICQlJ9Rw+Pj7MmDGDBw8ecPv2bVxdXenRo0eq+0ZERBAREUGbNm0A1UYjqVhsfHw8o0aN4vjx41hbWxMQEJDqMTK7nyXJVoISQlgBA6WUK3M4HouUtEAiuyv4UrKysuL111/n008/JSgoiOrVqz/1MTXNXLzbeScnIjFZICc9/Y3/UkqGDh3K9OnTn3gvPDycqKgo4uPjiYmJSW4fkVobCyklnTp14tdff031POm1nkhqO2FlZfVICworKysSEhKQUjJhwgTeeuutR74vJCTkkf2tra2Jjo5+4vgxMTGMHDkSPz8/KlasiLe3d5ZbZySZPXs2ZcuW5cSJExgMBhwcHJ5qP0uS7jUoIUQRIcQEIcQ8IURnY8280cAF4KXcCdH89u/fT+nSpZP/Qntar776KlZWVixevDhHjqdp+clzzz3HmjVruHnzJgC3b99ObkHx1ltv8dlnnzFo0KBHyhsdOnSI4OBgDAYDq1evplWrVjRr1ox///03uRXG/fv3c2zU8Pzzz/PTTz8lX9O6evVqcrxpSdneIikZlSpViqioqAxX7RUrVoxixYqxd+9egOS2IQB3796lXLlyWFlZsWLFChITE584X3r7WbKMRlArgDvAfuANVC09AfxPSnnctKFZjgMHDtCsWbMc64xboUIFXnjhBZYsWcLkyZOxtbXNkeNqmjlNajspR45Tt25dpk6dSufOnTEYDNja2jJ//nx8fX2xtbXl5ZdfJjExkRYtWrBz506srKxo3Lgxo0aNIjAwkPbt29O7d2+srKxYunQpAwcOTF6UNHXqVGrVqvXUMXbu3JmzZ88mTxk6Ozvz888/p9ut1svLi+HDh+Po6Mj+/fsZNmwY9erVw8XFJXk6Mz1LlizhtddeQwhB586dk18fOXIkL774IsuXL6dLly7JI0M3Nzesra1xd3fHy8srzf0sWbq1+IQQJ6WU9Y2PrYFQoJKUMntj0RyQ27X4bt++TcmSJZk2bRoTJkzIseNu3LiRXr16sWHDhlxpnaxp+dWuXbuYNWsWmzZtMnco2mNMXYsvuRWslDIRuGLO5GQOBw+qJsA5cf0ppa5du1K6dGmWL1+eo8fVNE3LLzJKUO5CiEghxD0hxD3ALcXzyNwI0Nz279+fPIWQk2xtbRk4cCB//vknd+7cydFja1pB0q5dOz16yqfSTVBSSmspZREpZWHjZpPieZHcCtKcDhw4gJubG87Ozjl+7CFDhhAXF8fvv/+e48fWNE3L6zJaxXdECPGNEKKLEMLy1yTmMIPBwMGDB3N8ei9Jo0aNqF27NitWrDDJ8TVN0/KyjKb4mgLrgHaArxBiixBijBDi6ZfB5AFnz54lMjIyR27QTY0QgiFDhrB3716Cg4NNcg5N07S8KqMpvgQp5S4p5UdSyqaopeb3gKlCiKNCiO9yJUoz2b9/P5DzCyRSGjRoEKBqj2mapmkPZalYrJTympTyJynlS4AnkK8rSRw7doyiRYtSs2ZNk52jcuXKtGvXjhUrVmRYLFLTNK0gyTBBCSEqCCE+FEJsFEIcFkLsNo6cuqJu4M23AgMDqVmzZo7doJuWIUOGcP78eQ4dOmTS82iapuUlGS2SWAL8BMQCXwADgZHA30AXYK8Qoo2pgzSXoKAgatSoYfLz9O3bFwcHB71YQtM0LYWMRlBfSSk7Sym/lVLuk1IGSilPSSnXSilHoxZP5F4HrlwUHx9PSEhIrhRzLVKkCL169WLVqlXExcWZ/Hyapml5QUaLJE5l8H6clDIwZ0OyDBcvXiQxMTFXRlCgpvnCw8PZtm1brpxP0zTN0qVbLFYIcRJI7cq9AKSU0s0kUVmAoKAggFxrh9G5c+fk0ke6Np+maVrG1cy750oUFiipRH9ujaCSqjR///333L59mxIlSuTKeTVN0yxVRlN8F9PbcitIcwgMDKRQoUK4uLjk2jmHDh1KXFwcq1evzrVzapqmWapM3QclhGhmXGIeJYSIE0Ik5vdisUndbk29xDwlDw8P3NzcWLp0aa6dU9M0zVJl9kbdeagl5ucBR1RFifmmCsoSBAYG5tr0XhIhBF5eXhw6dIgzZ87k6rk1TdMsTaYrSRhX61lLKROllEtQ90HlSwaDgQsXLuTaAomUBg0ahI2NDcuWLcv1c2uaplmSzCaoB0IIO+C4EGKGEOLdLHxvnnP16lViY2NzfQQFUKZMGbp27cqKFStISEjI9fNrmqZZiswmmSHGfUcB94GKwIumCsrccnsF3+O8vLwIDQ1lx44dZjm/pmmaJchUgjKu2ouRUkZKKSdLKd/LrzfoQu7fA/W47t27U7JkSb1YQtO0Ai2zq/haCiF2CCEChBAXkrYMvucnIcRNIUSq1SiEEO2EEHeFEMeN26fZ+QFMITAwEFtbWypWrGiW89vZ2fHyyy+zfv163Q5e07QCK7NTfIuBr4FWQOMUW3qWkvFCij1SSg/jNiWTsZhcUFAQVatWxdra2mwxeHl56XuiNE0r0DKboO5KKbdKKW9KKcOTtvS+QUq5G7j99CHmPnMsMX9cgwYNqF+/vp7m0zStwMpsgvIRQswUQjQXQjRM2nLg/M2FECeEEFuFEK5p7SSEeFMI4SeE8AsLC8uB06ZNSmkRCUoIwdChQzl48CBnz541ayyapmnmkNkE1RTVQXca8JVxm/WU5z4KVJZSugNzgfVp7SilXCil9JRSepYuXfopT5u+sLAwoqKizLZAIqVBgwZhbW3NkiVLzB1KgaObG2ua+WVULBYAKWX7nD6xlDIyxeMtQojvhBClpJS3cvpcWWHuJeYpubi40L17d5YuXcpnn32Gvb29uUPK82Ji4Px5CAh4uAUFQUQE3LuntshIiI8HOztwcAB7e/W1RAkoXx6eeUZt5ctDzZpQpw6ULQu5WBVL0wqETCUoIURRYBKQ1D3XF5gipbyb3RMLIVyAG1JKKYRoghrNpXtdKzckJShLGEEBjBgxgg0bNrB27VoGDhxo7nDynNhYOHgQfHxg5044cABS9oR85hmoUQOqV4fChdVWpIhKTrGxaouJgehoCA+Ha9fg2DG4cQMMhofHKVZMJSpXV2jcWG316oGtba7/yJqWb2QqQaHavp8CXjI+HwIsAfqk9Q1CiF9RHXdLCSGuoBKcLYCUcgHQFxghhEgAooEBUpp/YiUoKAgrKyuqVKli7lAA6NSpE9WqVeP777/XCSqTEhJg+3ZYsgS2bFHJRQho0ADeeQc8PaFWLTX6cXbO/jlCQ+HcOTh79uH2xx/w449qHwcHdc4WLaBtW2jdWiUyTdMyR2QmJwghjkspPTJ6LTd4enpKPz8/kx1/0KBB7Nu3j+DgYJOdI6tmzJjB+PHjOXnyJPXq1TN3OBbrzBmVlFasUCOcUqXgpZegUyeVIIoXN30MUsKFC3DoEBw+/PBrXJxKkh4e0K4dPP+8isnBwfQxaZq5CCGOSCk9s/39mUxQ+4EPpZR7jc9bArOklM2ze+LsMnWCatasGc7Ozvz9998mO0dWhYWFUaFCBd58803mzp1r7nAszvHj4O0NGzaAjQ106wZeXvDCC2qqztyio9U0o6+v2vbvV9OGjo7QoYOKs2tXqFrV3JFqWs7KrQTlASwDiqLavd8GvKSUJ7J74uwydYIqVaoUffv2ZcGCBSY7R3YMHjyYP//8k6tXr+Kc3XmpfObUKZWY/vgDihaF996D4cOhTBlzR5a+6GjYtQu2boXNm9WIC6B2bZWsXngBWrVSizM0LS972gSV2Vp8x43Lwd2A+lLKBuZITqYWERFBeHi4xSyQSGnEiBFERkby66+/mjsUs7t5E4YMATc3+Osv+PRTCAlRXy09OYEaOXXtCt9+C4GB6jrW7NlQsSLMmwcdO0LJkvDii/DLL3A320uRNC1vS3cEJYQYLKX8WQjxXmrvSym/NllkaTDlCOrIkSN4enqydu1aevfubZJzZJeUEnd3d2xsbDhy5Eiudvq1FFLCr7+qhQ737qkR04cfquXf+UVUlFpxuHmzmrK8fl1NU3bsqBJWz57q2pqm5QWmHkE5Gb8WTmXLd/NMlnQP1OOEEAwfPpxjx45x+PBhc4eT665cUb+cBw1Sq++OHYPp0/NXcgK1qrBHD1iwAK5ehb17YdQoOH0aXn8dXFzguedg/ny15F3T8rN0E5SU8gfjw7+NbTaSN+Af04eXu5LabFSrVs3MkaRu8ODBODk58f3335s7lFy1cqW6v2jnTjUVtncv1K1r7qhMz8oKWraEr76C4GA4cgTGj1eJa9QodaNwixYP39e0/CazpY5SWzqW75aTBQYGUq5cOZycnDLe2QyKFCnC4MGDWbVqFbdv58k6vFmSkADvvguDB4O7O5w8CWPHghmLzJuNENCwIXz+Ofz3nxpRffaZWnDxwQdQrZp6f+pUdT+WpuUH6SYoY3HY94HSQoj3UmzeQL77NWEJRWIzMnLkSGJiYli4cKG5QzGpsDDo3BnmzIExY+Cff9QvYU2pWxcmTlRTnUFBMHOmuqfq//5PvVe3rnp8/LiuK6jlXRmNoOxQ15psePT6UySqEkS+EhQUZJEr+FJyc3Ojc+fOzJkzh5iYGHOHYxJHj6pqD/v2wbJlKknpkkFpq1ZNjaL27VPX6ubOVbUBp01TlSxq1IBx49S9WDpZaXlJRtegfI3Xm5o9dg3qaynl+VyKMVfcv3+fa9euWfwICmDcuHHcuHGDFStWmDuUHLd5s7ruYjCoa02vvGLuiPKW8uXV9SkfH1WKaeFCtahk9mxo1gwqVVIj0j17IDHR3NFqWvoyew3qgbEf1BYhxM6kzaSR5bILxrslLX0EBdChQwcaNWrEzJkzScxHv2XWrIH//U8tiDhyRI2itOwrUwaGDYNt29S9Y8uWqetUP/wAbdqoZDZiBPz9t6rermmWJrMJaiXwH1AVmAyEAPlqrfPFixcBqJoH6s0IIRg/fjznz59n/fr15g4nRyxfDv37Q5Mm6npTXrjhNi8pXlyNRjdsUNf3Vq1SSWrFClWr0MVFlYdat07di6VpliCzCaqklHIxEG+c9nsN6GDCuHLdpUuXAKhUqZKZI8mcPn36UKNGDb788kssoAj8U/nhBxg6FNq3V5UhihY1d0T5W+HC6o+B335TyWrdOlVeaeNG6NNH3QjcrZv676LvtdLMKbMJKmkCIFQI0U0I0QDIV7dIXrp0CVtbW8qWLWvuUDLF2tqaDz74gMOHD+Pr62vucLJtzhxVP69bN9i0CSx0hX++5eioplWTKsD7+MDIkWop+/DhahqwSRO1fN3fXy+y0HJXZovFdgf2ABVR9z8VAbyllH+aNrwnmarU0csvv8zBgweTb9bNC6Kjo6lSpQoNGzZk69at5g4ny5YuhVdffVhzzhIqj2uKlKp9ycaNalrw4EH1eqVKqo5g166qEnvhwuaNU7NsuVIsFrgjpbwrpTwlpWwvpWyEqmieb1y6dCnPTO8lcXR0ZMyYMWzbto0TJ/JW7d4tW+CNN9T1D52cLI8QarHKhAmqC3FoKCxapJatr1ypRl0lS6pp2Rkz9OhKMw1dScIoLyYoUFXOnZ2dmTFjhrlDybSDB6FfP1Ud4o8/dHLKC1xc1B8U69dDeLiaCnz3Xbh9W5VfcneHChVUvcDff1f7aNrT0pUkgISEBK5evZonE1Tx4sUZPnw4q1at4mweqHFz7py63uTiokZReooo77GzU12Bv/wSTpxQNwcvXqzqAv7xh+piXLq0Gm299566thgZae6otbxIV5IArl27hsFgyJMJCmD8+PE4OTnx8ccfmzuUdIWGqlbnVlawfbuqdqDlfeXLw2uvqZHTrVvw778wZYpa2v7dd6o6e4kS6kbhCRNgxw548MDcUWt5QWYXSVSWUl4UQjgDSCnNdqeEKRZJ7N27l9atW7Nt2zaef/75HD12bvn888+ZOHEi//77Ly1atDB3OE+IjVX33Zw+rbrJ6ptwC4aYGNXifudOtR06pIoA29pC48aqakjLlmr0Vbq0uaPVclpuLZIoLIQ4BpwGTgshjggh6mX3pJYmr90DlZqxY8fi4uLC+PHjLfK+qNGj1S+n5ct1cipIHBzUQorPPlMjqzt3VKv7sWPVooo5c9SCizJloFYttapz0SK1gtBgMHPwmtnZZHK/hcB7UkofACFEO+NrlvenejZcvnwZgIoVK5o5kuxzcnJi0qRJjBgxgk2bNtGjRw9zh5Rs0SK1TZigbgTVCi5nZ+jSRW2gRlh+fip57dsHf/6pbj8ANUXYrJkaaXl6QqNG8MwzZgtdM4PMTvGdkFK6Z/RabjDFFN/bb7/NqlWrCM/jS4/i4+NxdXXF1tYWf39/rC2gcdLBg2pqr107tSjCAkLSLJiUEBCgktW//6ol7mfPPhxNlSunklXS1qiRvpZpyZ52ii+zI6gLQoj/A5LKZw8GLmT3pJYmry4xf5ytrS3Tpk2jX79+LF++nFdffdWs8dy4oW7CfeYZda+TTk5aRoSAZ59VW9I/3/v3VV8rPz+1HTmiVgYm/W1doYJa5l6//sPt2Wf17Qv5QWZHUMVRRWJbARJVVWKylPKOacN7kilGUO7u7lSpUoUNGzbk6HHNQUpJs2bNuHbtGgEBATg6OpoljoQE6NhRjaD27VNLjjUtp9y7p5o1JiWskyfVSCshQb1vawu1az+atOrWVZUw9B9KucfkIyghhDWwVkrZPrsnsXSXLl2iTZs25g4jRwgh+PLLL2nfvj3ffPMNH330kVni+Owz8PVViyJ0ctJyWuHCauo45f9t4+LUfXYnTz7c9uxRo/ck9vaqP1bt2g9HakmbLlJseTJMUFLKRCGEQQhRVEp5NzeCyk2RkZFERETkiym+JO3ataNHjx589tln9O/fP9dbiPz7ryouOnQoDBmSq6fWCjA7u4ejpZQiIuDUKTXCOndObSdOqCruKduplS2rVhJWqwZVq6qvSY9dXNT9e1ruyuw1qCjgpBBiB3A/6UUp5TsmiSoX5YcVfKmZN28erq6uDB8+nG3btiGEyJXz3r0LgwZBlSrw7be5ckpNS1exYtCqldpSiouDCxceJq1z5+D8edWP7OrVR2sLOjiof9NJCatSJXXtq2JFtT3zjL7mZQqZTVBrjVu+kx/ugUpNpUqV+PzzzxkzZgwrV65k8ODBuXLet99WpW/27oUiRXLllJqWLXZ2aqqvdu0n34uJgYsXIThYJbGUX//9V/0hlpIQagRWseLDxPXMM+q1pM3FRd2MbGubOz9ffpCpBCWlXGbqQMwlvyYoUMvnf/nlF9599126dOlCqVKlTHq+lSvVNmWKun9F0/IqB4eH16ZSc++e+kPs8mW1pXx87hz8/bfaJzUlSz5MWCkTWNJrpUqp0lAlSqjrYgV5ajHdBCWE+BN1Q+42KWX8Y+9VA7yAECnlTyaL0MQuXbqEtbU15cqVM3coOc7a2ppFixbRsGFD3nvvPZYvX26yc4WEqEZ3LVuqG3I1LT8rXBjq1FFbWu7fV7da3LgB168/fJzy+aFD6vH9+6kfQwh1w3JSwkpvK1pUzVoUKaLiK1Ik74/WMhpBDQPeA+YIIW4DYYADUBUIBOZJKVNdmy2E+AnoDtyUUj5RFkmoiyLfAC8ADwAvKeXR7P4g2XXp0iUqVKhgETe1mkL9+vUZP348n3/+OYMHD6Zz5845fg6DAV55RT3++WewyezEsablY05ODxdaZCRlMrt1S7UxeXy7c0e1MTl/Xj2PiMi4B5eDw6MJK7XHzs5QqFDWNltblTxNLVP3QQEIIaoA5YBoIEBKmW49YiFEG9TiiuVpJKgXgNGoBNUU+EZK2TSjOHL6Pqi2bdsipWT37t05dkxLExMTg7u7O/Hx8Zw8eRKnHO6rPncuvPOOKlEzdGiOHlrTtDQkJqprYbdvq8R1966aVoyMVFvKx48/T/k4Njbr57a2fjJpOTqqZfz29ur6nr09bNiQO5UkkFKGACFZ2H+3MamlpRcqeUnggBCimBCinJQyNLPnyAmXLl2iZcuWuXnKXOfg4MDChQtp3749o0eP5qefcm5GNiRETel16fJwFKVpeZH3Lm+823mbO4xMs7Z+OL1Xo0b2j5OQANHRqgVKRlt6+92/r1ZGxsaqBJidxPc4c07GlAcup3h+xfhariWoxMRErly5ki8XSDyubdu2fPLJJ0ydOpXWrVvnSBkkKWHYMDXU/+GH3Bnya1pGsptoJvtOzlMJKj2xcbFEPojk3oN7REVHce/BPe7H3Od+zH2i46KJjosmJj7miS02IfbhlhhLXEIccYlqSzAkEJ8YT6IhkXgZT4J9Aom2iSQUSSDRkEiiVJtBGjDIRKRMBP+n+znyxNUCIcSbwJuQs6vtrl+/TkJCQoFIUADe3t7s27ePkSNH0rBhQ9zdn67W75IlarXS99+r+0I0zRJYYqIxGAxEREVwI+IGtyJvERYZxu17t7l9/7ZKJNH3uBd7j6iYKO7H3udB3AMexD8gOiGamMQY4gxxxMk4EkggnngSSCBRJJJopTZpLZFWasOazDdSyiKB8fACrAXYCPXcRoCNfLg5GMA6B9qlZCpBCSH6AJullDkwaEt2FUh5d2wF42tPkFIuRK0mxNPTM8eaHeXnJeapsba25pdffqFBgwb069cPPz8/imTzZqVr11Q777Zt4c03czhQLd+zlOk0713eTPadnPxcTFbTAJPaTsK7nTdxcXFcD7/OlVtXuHb7GtcjrnMj8ga3om4Rfj+cO9F3uBt7l6j4KB4kPCA6MZpYGUsccSRYJ5BonYjB1gC2ZC1pWAE2qnSZsBJYJVphbbDG1mCFvcEaJ4M1jgZbChmscEq0wilR4JQgKJQIhRIMOCVICsUZcIxPxDHOQKG4BBxiE3CMScA+QWKfCHaZ2GwNAjs7R2ztHLB2KISVYyGEg4O64JS0pfNc8HRdvjNbLHYJ0AHYDaxGLTtPyMT3VQE2pbFIohswioeLJL6VUjbJ6Jg5uUhi9erVDBgwgJMnT1KvXr7pv5ih3bt306FDB3r37s1vv/2W5SoTUqomczt2gL//081/awWTmCyQk3KusebjiSZJUqJJTEzk1u1bBF0LIvhmMJduXeJqxFVuRN7g5v2b3Im9w8lCJylxpwQxxBBnFUeibSLS3jgiSU88WMVbYZ1gjY20wU7a4WDlgL2wx9HakULWhXCydcLZzhlnO2eKOBShqEMRSloXopTBhhIJUCLWQPHoBIo/iKNIZDRF7tzD9tZtCAtTy/ru3FGrIDL6fV2kiCqdUayYWndeuPDDpXqpbam95+T0cNXDU65Tz5V2G1LKV4UQtkBXYCAwXwixQ0r5RjqB/Qq0A0oJIa4Ak1B/RyClXABsQSWnQNQy81zvDVHQRlBJ2rRpw7Rp0xg/fjxz587lnXeyVrHqt99g40aYNUsnJ818DAYD129e5/TF05S6U4qxLmO5eOci62LXUSOyBnfj7zJz7UymbphKon0ipFXY3xqElfojLU7E4YQTZazKUNi6MMWsi1HMsRglC5WklHMpShcuTdliZSlXvBzlS5anQqkKFHVKUWU2Pl6tFb92TdVLunZNbZeuQWgo3AxSSScsLO1VBA4OquRE0larlloJkZR40tqKFs13pdozvcwcwJikuqCSSRsppWlLE6QiJ0dQo0ePZsWKFUREROTI8fISg8FA79692bx5M+vWrct0B967d1VpmAoVVDO5fPb/B82EMhrlpHT//n38g/w5Hnycc1fPEXwrmCuRV7gZfZNIQyQPrB+Q4JgATqgLI4+xfmCNXYIdTjhRxKYIJexLUMqpFGWdy1KuaDkqlqhI5dKVqeZSjWou1bC3tc942lFKuHlT1UBK2i5dUl8vX1aJ6ObNJ0c51taqRES5cqq3fZkyjyagxzcnp3yz4uhpR1CZneLrCvRHjYh2Ab8Bf2Vmmi+n5WSC6tWrF8HBwfj7P+VSkzzq3r17dOjQgVOnTrF9+/ZMtRx55x2YP1/dAd+oUS4EqZlc0i/mtH5B5/T1ooSEBGw/t+W3Zr9x8tJJzt04R0hECNdjrnPHcIf7dvcxOBuv26QkwSbWhkKJhShqVZSS9iVxKeRCxWIVqVa6GrXK1cK1kisrz61kSocp2QvuwQMICoLAQHVHbGCgKsCXlIweH/UULgyVK6tVQuXLqwJ8j2+lSxfYv+RyK0H9irr2tDWHF0pkWU4mqAYNGlC+fHk2bdqUI8fLi27dukWrVq0IDQ3F19cXDw+PNPc9ehQaN1YljebOzb0YNdNKuh6U1nWh7FwvunPnDgGBARwKOMSxi8c4F3aOy/cvE24I54HjA0hl7sUm1ganBCdKWpekXKFyVC1elWddnqVuxbq4V3OnUolK2FrnQO2e+HiVeM6eVUkoKREFBqppuZRKlVKlICpXfpiIkh5Xrqym1rQ0mfQalBCiBlBWSjnwsddbAtellEHZPbEluHz5Ms2bNzd3GGZVqlQpduzYQcuWLXn++efZu3cvNWvWfGK/xEQYPlz9MfjZZ2YIVLM4UVFRnD9/Hr8zfhwIOsDpG6cJiQoh3CqchCIJUIyHK9eKgZWzFUUSilDJrhLxtvF0r9SdehXq0aBaA2q71MbJLmcrnBAXp5LO6dNw5szDrwEBKkklKVtWXUzt2FF1M6xR4+GmuxiaVUaLJOYAqZX+jDS+l7kLFxbo/v37hIeHF7gFEqmpWLEif/31F61bt6ZTp078+++/lC9f/pF9Fi2Cw4dVtXL9R2Pel9by6qSvbSu3xfei7xPv13hQA67CtfhrPCj0AEqjrgMBlAGrklaUNJSkvGN5apaoiVsFN5rWbEr9CvUp61TWdH3JwsNVD/ijR9VXf3+ViJJ6wAuhRkKurtCjh/pap45agFC4sGli0p5aulN8QojDUsrGabx3UkpZP7X3TCmnpvj+++8/6tSpw8qVK3n55ZdzILK8z8/Pj/bt21O6dGm2bdtGrVq1ALUoqXZtaNhQ3ZibT67fFijpXUdKOcXn09aH4yePs/DKQhKDE7kYe5HYhrGqTWmKAY59oj0uNi7ULFaTBhUa0PLZlniU96Bi0YpYCRP2h5BSTcOlTEbHjqnrQ0kqVQIPD5WEXF2hbl31D9gxrWV8mqmYepl5sXTey9P/tQvqEvP0eHp68s8//9C9e3datGjBxo0badGiBePGqTpb8+fr5GRpMruAIWV1hbt37/LuhndpGt2UY6eOQSlw7uQMraD96vZqVFQIcAVrqS7ud63WlZY1WtKyRkvqlq5L6UKlc6dL84MH4OcH+/er7cAB9RcTqH+MtWqpHi+jRkGDBmorWdL0cWm5IqME5SeEGCalXJTyRSHEG8AR04VlejpBpa5Jkybs27ePrl278txzzzFx4naWL2/DJ5+k3nlUM6+MyvqEhoZy7NgxAPr07cOhi4e4Kq9CD1hyfQmUUfvdb3UfO+xwq+1G08pNmX9iPqdHnqZWyVpM3T01d6o+SKmqDyclo/374cSJh9N0NWpA587QpIkazru5qRtLtXwrowQ1FlgnhBjEw4TkCdgBvU0Yl8ldunQJKysrnnnmGXOHYnFq1KjBvn376N69FxMnOlOs2D0+/ljP0+emrC7tNhgMBAcHc+zYMY4dO8bRo0fxdfQl2j06eZ919ddBikn5Ng3b0KZ6G5qUb0Lj8o1Z4LeAyb6T8YtQU+iu37kC6j4lk5BSLWLw8YGdO8HXV3XvA1XJoEkT+PBDaN5ctWguXdo0cWgWK7PLzNsDSbWATkspd5o0qnTk1DUoLy8vdu7cmTyS0p70/fexjBxpDwxg6FAHvv3222zX7tPS5r3LW31NkZDSW9o9yWcSU3ancp/PIeAG0ARsnW2Jd1Ir1RysHIgxxNCsfDMOXD3w5PFSuVE2p0sRJbt4USUjHx+1XbmiXi9XDtq3V9N1zZtD/fq682U+kFuljnwAn+yexBJdunRJT++l4+5d8Pa2p0ULSYcOtZg27XN2797NypUrC/zS/JyWtJoutRGTlJIrV65w+PBhDh06xOHDh/Hz81OtQMsBw6DYzWJEl4omtsnDWxR71ulJq0qtaFWpFe5l3bGbasf+N/Ynv2+yBPS4iAj46y+17dypbnoFdX9R+/Zq69BBXUvSFzi1xxTYP1EuXbpE48apLlDUgKlTVbmwLVsEjRpNoUuX5xk8eDCtW7dm4sSJTJw4ERv9F26OSmvpN7vUZl3WmmdaPUOJkSWId4wnWqrpu7J1yiYno1aVWlFzbk3WvLTmkWNndZou29N6Uqr7jTZvhi1b4N9/1U10xYqp0vdjx6qk5OoKViZc7aflC1mqxWcJcmKKLzExEQcHBz788EOmTZuWQ5HlH+fPq98fQ4bA4sUPX797925y/cImTZrwzTff0KxZM/MFmoelVZcOUAmpHVTfUB2X5i7EVYgjRIQQFhsGQNViVelYrSPPVX0Ov2t+zOw8M0t17pLOn2MLH+7fV6OjLVvUljRt7uEBL7ygtqZN9ZRdAZQrpY4sSU4kqIsXL1KlShUWLVrEG2+kWZC9wOrZE3btUvc5urg8+f7q1asZO3Ys169fp3///kyfPp2qVavmepx5SVxcHP7+/hw+fDh5O3PmDIY2BlXhEuh9sjelG5UmokwEv135Lfl7SxUqxXNVn1NbteeoVrxauufKlem78HDYsAHWrFHJKTZWFTnt1Am6dYOuXVVtOq1Ay5VrUPnNhQsXAPQv1VRs3w5//glffpl6cgLo378/3bp1Y8aMGcyaNYt169YxZswYJkyYQPHixXM3YAuUmJjI2bNnk68XHT58mBMnThAXFweo8lKNGzemT58+TLF6uNhhp+dO7sbexfaaLVWLVeXtxm/TsVpH6petb9qbXzPrxg1Yv14lJR8fNXVXtSqMGKGSUuvWYG9v7ii1fKRAjqB++uknXn/9dS5cuKCTVArx8eDurkqYnT6dud81V65cYeLEiSxfvhxHR0eGDh3K6NGjqVOnjukDtgBSSoKCgh5JRkePHuX+/fsAFC5cGE9PT2JbxDLWfSyenp5EOESw+fxmNgVs4uDVgwA42ToxoN4AutXsRsdqHSlsn/1l/Tk6fXftGqxdC3/8Abt3g8Gg6tX17au2Bg304gYtTXoElQ0XLlzA2tqaihUrZrxzAbJggSrwvGFD5v8QrlChAkuXLuW9995j9uzZLF68mO+//57OnTszevRounTpkqcXU6T8ZR8fH8+5c+f4dOen1LhSg+PHj+Pn58edO3cAcHBwoEGDBriOcmV0vdF4enpSq1YthBBYTbGiRdEWjNs4jpCIkCfOcz/+PhWKVKB3nae/vfCpk9Pt26or5c8/w759auFD3bowcaJKSvXq6aSk5YoCOYJ6+eWXOXDgQPJUn6a6StesCZ6eakVwdn//3Lx5k4ULF/Ldd98RGhpKiRIl6NmzJ3369KFTp044ODjkbOAmcufOHU6cOEF73/a8evFVTpw4walTp9Q0nTfYT7fH1dUVT09PGjduTOPGjalbty62traIyYLETxPZd3kfa86sYe3ZtVyOvIytlS0dq3XkxTov0r1Wd8o6l8295d4ZiYtTCxxWrIBNm9RzV1cYMABefFEVVtW0LNIjqGy4cOEC1aqlf6G5oJk0Ce7dg9mzn+6P4zJlyjBx4kTGjRvH5s2bWbt2LevWrWPp0qU4OzvTuXNnWrZsSYsWLWjQoAH2ZrxmERcXR3BwMAEBAQQEBHDu3Lnkx6GhoWonb9i8eTMeHh6MGTMGd3d3BgcOJioq6omRYaIhkV0huwCo8HUFQqNCH3k/3hDP1sCtNCnfhLLOZXPhJ8yAlKrz5IoVsGqVWvhQpgy8/bZawunhoUdKmlkVyBFU2bJl6dWrFwsXLsyhqPK2kyfV7yJTNSKMi4vDx8eHtWvX8tdffxESEgKAvb09jRo1wt3dnZo1ayZvVatWxc7O7qnP++DBA27cuMGVK1ce2YKCgggICODChQskJiYm71+qVCmeffZZ7je+z/FixzN1jraV2/JNl28YsXkE+6/sf+L9j1p+xPSO09McKeV0t9pMuXwZli1TiSkgABwc4H//U0mpc2e9HFzLMXqZeRZFRUVRuHBhpk2bxoQJqbW6KlikVH3ajh1T9z/lRiHo0NBQ9u/fz759+9i/fz9nzpwhIiIi+X0hBMWLF6dkyZLJW9GiRbGxsXlki4+PJyYmhujoaKKjo3nw4AG3b9/m1q1bhIeHExMT88S5CxcuTPXq1alVq9YjW82aNSlRosQT+6dMLCnbUlx97yor/Vcy7u9xANhY2fBCzRcYVH8Q/df0fyIZmX0qLzERtm2DH35QN9EaDOrG2SFD1HUl3ZhPMwE9xZdFwcZSK3qKT9mwQd3GMndu7nUpKFeuHH369KFPnz6AWgkXHh7O+fPnOX/+PBcuXCAsLIzbt28THh5OaGgo586dIyEhIXmLj4/H1tYWBwcHHB0dk7dKlSrRoEEDSpUqRcmSJSldujQVK1akQoUKVKhQgSJFijzVqGXFiRUAVJxdEYM0ADCv6zz61+tPqUKqj3n/Nf2f+D6TFVzNyLVr8NNPquPkpUuqe+xHH8Ebb6gl4ppmwQrcCGrjxo306tWLQ4cOFfhSR7GxanGWoyMcP15wZnayMpppu6Qtuy/tztS+SVUbzDJtl5LBoDpLLlgAGzeq0VPHjvDWW9CrF9jami82rUDRI6gs0jfpPjRnDly4oFbtFZTklFl3ou/ws//PhEeHA1DEvgiRsZHseXUPLSu2TG7Wl1qyM1tyunNHjZR++EH9hy1VCt57D958U/VS0rQ8xgJuT89dFy5coHDhwpQs4F03Q0NVQdiePVV1mvzOe5c3YrJILsCa9Dip1QWoqcb9l/fjtd6L8l+X551t71DIthA/9viRa+9dA6BVpVa500k2K86dUyvvKlSA8ePV119+Ua0sZszQyUnLswrc381JS8wt7pdMLvv4YzXF99VX5o4kd3i3ezjt9vioJyImgp/9f2bhkYWcvHkSZztnXnF/hbcavUWDcg2S90vtOpLZri1JCf/8o4bBmzeDnR0MGqSqhbu5mScmTcthBS5BBQcH8+yzz5o7DLM6fBiWLoVx4wr2H9cB4QHMOTCHZSeW8SD+AQ3LNeSH7j8wsN7AVEsNpVoVPLen82Ji1Ohozhx1f0CZMuDtDcOHqwUQmpaPFKgEJaXkwoULdOnSxdyhmI2UMGaM+l32ySfmjib3SSl5xe0Vevzag00Bm7CztmNQ/UGMbDwSz2eyfS3X9MLCYN48+P579djNTa3OGzhQ3cekaflQgUpQ169fJyYmpkAvMf/1V9i/X/V5Kkjd22MTYll1ahWzD8zmxI0TlC5UmkltJzHCc4RlVHVIy+XLMGuWWvwQHQ3du8O776qmfwV8mlrL/wpUgiro90Ddv6+m9Ro1Ai8vc0eTO249uMUCvwXMPzyf61HXcS3tyo89fmSQ2yAcbCx45PHff6rnyc8/q+eDBqkFELomnlaAFKgEVdCXmM+YAVevqrJrebHbdlbuL7oYcZEZ/87gp+M/EZMQQ5caXXiv2Xt0rNbRshfIHDkC06erFhcODqrX0vvvQ+XK5o5M03JdgUxQVapUMW8gZnDxokpQAwZAq1bmjiZ7JvtOzjBBBd4OZPqe6Sz3X45AMNR9KO82f5e6pevmTpDZIaVqYTx9OuzYocoOTZigLhaWKWPu6DTNbApUggoODqZ8+fJ5puVDTnr/fXXJ4ssvzR2JaZwJO8O0PdP49dSv2FnbMcJzBB+2+JCKRS2451fSUnFvb/j3X5WMvvhCrcjTtfE0zbQ36gohugghzgkhAoUQH6XyvpcQIkwIcdy4vWHKeApqB91//lENUT/+GCpVMnc0WZPRDbbHQo/R97e+uH7nyvr/1vN+8/cJHhPMt12/tdzkJKUqRdS6tbpLOiRErdALCVHXmXRy0jTAhLX4hBDWQADQCbgCHAYGSinPpNjHC/CUUo7K7HGfphZfxYoV6dChA8uWLcvW9+dFSW3cY2NVG/e8PHhMeYPtwSsH+Wz3Z2w+v5mi9kV5p+k7jGk6hpKFLLhCiJSqMq+3N+zdC+XLq6m8N97IfAtjTctDLLkWXxMgUEp5AUAIsQroBZxJ97tMJDY2lqtXrxa4FXzz5j1s456Xk1OSkzdOMtFnIhvPbaSkY0mmtp/K203epphDMXOHljYpwcdHJaY9e+CZZ9R/mNdfzx//UTTNREw5xVceuJzi+RXja497UQjhL4RYI4RIdU5GCPGmEMJPCOEXFhaWrWAuXryIlLJATfHduKF+J3bpAj16mDuatLVb2i7DfYJuB1G/TH3cF7jjG+LL1PZTCRkbwidtPrHs5LRrl+q79NxzEBSk+poEBanaeTo5aVq6zL3Y+E+gipTSDdgBpDr3JqVcKKX0lFJ6li5dOlsnSlrBV5BGUB99pO7t/OYby76n0/eib5rvXbt3jRGbRlB7fm0CbwcyruU4Loy5wCdtPsHZzjkXo8yiw4fV9aX27R9NTKNG6cSkaZlkyim+q0DKEVEF42vJpJThKZ7+CMwwVTAFLUEdOPCw3l6tWuaOJuvCH4Tz5b9fMvfQXBIMCbzZ8E0mtplIucLlzB1a+s6ehYkT1X1MJUuqarwjRqimW5qmZYkpE9RhoKYQoioqMQ0AXk65gxCinJQy1Pi0J3DWVMEEBwdjb2+Pi4uLqU5hMRITYfRoKFdO/a60RO2Wtntk5JS0Sq9VxVY8X+N5Zu6byb3Yewx2G4x3O2+qFbfwPyxCQmDyZFi+HJyc1Nzqu+8WrHpSmpbDTJagpJQJQohRwHbAGvhJSnlaCDEF8JNSbgTeEUL0BBKA24CXqeJJWmJulRdLKGTR99+Dn58qel34yaLcFmGX167kx2KyIPHTRJafWM4nOz/h/3z+j//V/h9T20/FtYyr+YLMjBs34PPPVfdaKyvV7mLCBNUsUNO0p2LSG3WllFuALY+99mmKxxOACaaMIUlSH6j87upVdb9T586qakRe4bnQk2PXj9G0fFN+7/c7LSq2MHdI6YuIUEVc58xRLTBeew3+7/+gooXee6VpeVCBqCSR1GajVV6t8ZMFY8eqe5+++86yF0aA6sc0bsc4AMKjw/n1xV/p79rfsmvlPXigloh/8YVqsd6/P0yZkjcv9GmahSsQCerOnTtERkbm+yXmmzfDmjVqxql6dXNHk7bwB+FM8Z3Cd37f4WjjyPTnpjOm6RgcbS14IUF8PPz4I3z2GYSGQteu6oNu0CDj79U0LVsKRIIqCCv47t9Xt9bUrQsffGDuaFIXlxjH/EPzmbJ7CpGxkQxrOIzJ7SZbdj8mg0E10fr0U7hwQVXaXb1alSnSNM2kdILKJyZPVhXL9+wBOztzR/Okree3MmbbGM7fPk+XGl2Y2Wkm9crUM3dYaZMSNm1SbYdPnlT1ojZvViMnS56C1LR8pEAkqKRGhfl1is/fH77+WpV0s7TLbIG3A3l3+7tsCthErZK12PLyFrrW7GrusNLn66tW4u3fDzVrqgZa/frlzSZampaHFYgEdfr0aVxcXChsqWuun0J8PLz6KpQoYVmtNKLiovh89+d8feBr7KztmNlpJu80fQc7awsc3iXx81Mjpr/+UoVcFy5UrYdtbc0dmaYVSAUiQR04cICmTZuaOwyTmDoVjh5VhQtKlDB3NGrF5C8nf2Hc3+O4du8aQ92HMv256ZZdAeLsWbVE/I8/HlZ/GDlSlyTSNDPL9wnq9u3bnD9/nldffdXcoeQ4Pz+1kGzIEOjd29zRwNHQo4zeOpp9l/fh+Ywnf7z0B80qNDN3WGm7eFFVfFi+HAoVgkmT4L33dPUHTbMQ+T5BHTp0CCDfjaCio+GVV8DFBb791ryx3Hpwi0/++YRFRxdRqlApFvdcjJeHF1bCQq/Z6OoPmpYn5PsEdfDgQYQQNG7c2Nyh5KiJE9XM1PbtUKyYeWJIMCSwwG8B/+fzf9yLvcc7Td/Bu5235ba/iIiAmTNV9YfYWFX94dNPoUIFc0emaVoq8n2COnDgAK6urvlqgcTu3TB7tiqS3bmzeWLwDfFl9NbRnLx5kueqPse3Xb+lbum65gkmIw8eqGHml1+qJDVggKr+ULOmuSPTNC0dFjoHkzOklBw6dChfTe/du6cWllWtCjNM1pwkbVcirzBgzQDaLWvH3di7rOm3hh1DdlhmcoqLUzWfqldXU3gtWsCxY+rGW52cNM3i5esRVGBgILdv3843CUpK1SX84kV1q45zLvbri0mI4ev9X/P5ns8xSAOT2k5iXMtxFLItlHtBZFZcnGqG9fnncOmSqvrw+++Wd5OYpmnpytcJ6uDBgwA0a2bBK8myYM4c9Xv2yy9z73etlJJNAZt4d/u7BN0Jok+dPnzV+SuqFKuSOwFkRVwcLFkC06apxNS0KfzwAzz/vK7+oGl5UL5PUM7OztSta4HTT1m0Zw98+CH873/qa24ICA9g7LaxbA3cSp1Sdfhr8F90qt4pd06eFY8npmbN1E22nTvrxKRpeVi+TlAHDhzA09MTa2trc4fyVK5fV10dqlZVM1em/p17L/YeU3dPZfaB2TjYOPBV568Y3WQ0ttYWVlEhNvZhYrp8WSWmRYugUyedmDQtH8i3CSomJoYTJ07w/vvvmzuUp5KQoBadRUTAtm1QtKjpzmWQBlUFYsc4QqNC8fLwYvpz03FxdjHdSbMjOlolpi++UImpeXPVCkMnJk3LV/Jtgjp27Bjx8fF5foHERx+pBRErVoCbm+nO4xPswwc7PuBo6FE8n/Fkbf+1llcFIjwc5s9XDQPDwlRiWrwYOnbUiUnT8qF8u8w8aYFEXk5Q3377sCzc4MHZP473Lu803zt98zTdf+lOh+UduPXgFit6r+DgGwctKzkFB8Po0aqd+qRJavGDry/8+68eNWlaPpZvE9SBAweoWLEi5cpZcJHSdKxcCWPGqBp733zzdMea7Dv5iddC74UybOMw3Ba4sffSXr7s+CXnRp1jsNtgyylRdOSImt+sUUOtxhswAE6fhj//hDZtdGLStHwu307xHTx4MM8uL9+8Wd2M2749/PIL2GTyv5L3Lm+823mnu09UXBQz/53JrP2ziE+M550m7zCxzURKFir51HHnCINBtbuYORN27lSFWz/4AN55R7XA0DStwMiXCermzZuEhIQwatQoc4eSZf/+q3rjubnB+vVZ6/gw2XdycoLy3uX9yMhJTFajDSdbJ+7H3+cl15eY1mEa1UtUz8Hon8Ldu7BsmbrGFBCgktHMmfDmm7q6uKYVUBYyl5Oz8ur1p+PHoXt3dall69b0fy+nd10JwLudN3KSJG5iHEDyjbUNyjXgwOsHWN13tWUkp5MnVVHB8uXVnGaJEmpFyIULauSkk5OmFVj5cgR18OBBbGxsaNiwoblDybS//oK+fdUy8r/+gjJl0t8/abSU1khpYuuJVClWhal7pgJQqlAp5nWdxws1X0CY+9pNZKRqo754MRw6pIaJAwfC229Do0bmjU3TNIuRLxPUgQMHcHNzo1AhC6wTl4off4Thw8HVVV1/ykr3B+92D687icmCexPusfjoYmYfmM3FuxfxfMaT5hWas7LPSvMmJoNBlcNYuhR++01VGK9XT5VlHzJEdbLVNE1LId8lKIPBwOHDhxk0aJC5Q8mQwaA6jU+bpsrF/fbbkzNa7Za2Y5fXLiDt60qT2k5iuOdwACrNrsSdmDu0rNiSeS/Mo1vNbuZLTFKCv79akvjrr3DliqpwO2iQqnrbpIleiadpWpryXYLy8/MjMjLS4q8/RUWp6/+//qq+zpsHtqlUEvK96Jvhsdb/t54v9n4BQLsq7fiwxYc0r9g8p0POnKSktG4drFmjloXb2KgMPGMG9OwJTk7miU3TtDxFSCnNHUOWeHp6Sj8/v1TfS0hIoHnz5ly6dImzZ89SokSJXI4uc7Zvh7feUnVNp0+HcePSHkiIyQI56dH/Rvdi71HkiyK4l3XnxI0TONs5M8RtCGObjaVWyVq58BM8JiEBDhxQyw7XrVMLHIRQJdcHDICXXtLt1DWtABJCHJFSemb3+/PVCGr27Nn4+fnx22+/WWRyunUL3nsPVlz2prajN3v2QMuWT+7Xbmm7R0ZOSVN5pQuVJuxBWPLrJ26cAGBU41FM7zjdtME/7soVlWm3bYMdO9QycVtbVXboo4/USKls2dyNSdO0fCXfjKACAgJwd3ena9eu/PHHH+ZfqcbDG2fj49VlmHHj4M4dSJgoiPlIYm+f9veKyYLYibHYT7XnrUZv8fuZ37kdfZviDsV5sc6L/HjsRwyfGnLv57x8WfWa371blRk6d069Xr48dOmitk6dTFvNVtO0PEWPoFALI15//XUcHByYP3++SX5pZ6ZKw+Mm+07G4YA38+bB1avQuDH88w+4rSXN5HQl8gpbz28FoOQMtbJthf8Kej3bi4H1BvJ8jeexs7bjx2M/mi453bunygz5+cHhw2opeEiIeq9IEdWh9vXXVVKqV08vdNA0zSRMmqCEEF2AbwBr4Ecp5RePvW8PLAcaAeFAfyllSFbP891337F3716WLFmSZu291BJMWkkntddTVmlIT0KC+n2+YgXgAhMmwHPPQZtPvfk1dDJua9V+SdN2wxsNx8PFg/1X9vPnuT+5HXM7+VhRcVEAjGk6hmnPTXvkPJPaTsowlgzFxEBgIJw6pRYzJG3nz6vFDgCVK4OnJ4wdq+rfublBHu+vpWla3mCyKT4hhDUQAHQCrgCHgYFSyjMp9hkJuEkphwshBgC9pZT90zvu41N8ISEh1KtXj1atWrF169Y0RxWpLTZI7bWs7puQoH6f79ypLsVsjfEmrvmTxVlHNxlN1xpdOR12mg93fEjT8k05HXY6OQmVcSpD8wrNaVWpFV1rdKXe9/VSPV+mSanaU1y79nC7elUtYLhwAYKC1PMkVlaqKKurK3h4qOFeo0YZ3zGsaZqWBkue4msCBEopLwAIIVYBvYAzKfbpBXgbH68B5gkhhEwna0ZFRrF80TJCQoK5ePESR/wOU8yxMB+9N56zx08jAaREIpEGA9L4GOD4IT/1XEqkNABweN9+pFT7ISXxCeo8O7bsYunlJfxyfXnyuZNGPQ3vDaTs+T5cC4vm5p1YEqwfYHCIpIhLBA1rRVLEvh/WxW6z9cY/lLYvQVjsbeYemsvcQ3OTj+UUL/Cq1IsmJevToqQH1ZwqquSamAiXY9ROPj6qa2xcnPqatMXEqEUJERHqolZExJOPIyLUsR73zDNQvbpazFCt2sOk9OyzWSv8p2maZmKmHEH1BbpIKd8wPh8CNJVSjkqxzynjPleMz4OM+9xK87jPCMlbJgk5R9gnQIloKPUATpaFN/2gYiRUiISa4VDnFnzbFLx3pX8c73YZ74OjIxQr9nArXvzR52XKqIRUvrz66uKik5CmabnGkkdQOUYI8SbwJkDRIoXoE9AaR3sHbGxsEAgQAgE8/F+BTHok1PNZ5VbywbXBWImHe35Zbjnjr3shZNJ3CaysBVPLLMb79jCsrcDGSmDvAO/ZLeRH++HY2wmc7W0oZGWHo7DD0cqOQsKOYtZOFLdywtH64eoH7zvr8O7f5+EiAuNXb4A3H30t+auVFdjb421vD5Ps1WqK1LaiRdNeaaFpmpYPmHIE1RzwllI+b3w+AUBKOT3FPtuN++wXQtgA14HS6U3xpXejbrrxPOU1qOys4tM0TSvInnYEZcp2G4eBmkKIqkIIO2AAsPGxfTYCQ42P+wI700tOTyO1VW9prYRL7XWdnDRN03KXSW/UFUK8AMxBLTP/SUr5uRBiCuAnpdwohHAAVgANgNvAgKRFFWnJ7ghK0zRNy10WfQ1KSrkF2PLYa5+meBwD9DNlDJqmaVrelC876mqapml5n05QmqZpmkXSCUrTNE2zSHmumrkQIgy4aO44UlEKSPMGYwuWV+MGHbs55NW4QcduDs9KKQtn95vzxI26KUkpS5s7htQIIfyeZrWKueTVuEHHbg55NW7QsZuDEOKpllzrKT5N0zTNIukEpWmaplkknaByzkJzB5BNeTVu0LGbQ16NG3Ts5vBUcee5RRKapmlawaBHUJqmaZpF0glK0zRNs0g6QWWDECJECHFSCHE8aRmlEKKEEGKHEOK88Wtxc8cJIIT4SQhx09gcMum1VGMVyrdCiEAhhL8QoqH5Ik8zdm8hxFXjZ3/cWJA46b0JxtjPCSGeN0/UIISoKITwEUKcEUKcFkKMMb5u0Z97OnHnhc/cQQhxSAhxwhj7ZOPrVYUQB40xrjZ2VkAIYW98Hmh8v4oFxr5UCBGc4nP3ML5uEf9eUsRvLYQ4JoTYZHyec5/5wxboesvsBoQApR57bQbwkfHxR8CX5o7TGEsboCFwKqNYgReArYAAmgEHLTB2b+CDVPatC5wA7IGqQBBgbaa4ywENjY8LAwHG+Cz6c08n7rzwmQvA2fjYFjho/Cx/Q3VJAFgAjDA+HgksMD4eAKw2R9wZxL4U6JvK/hbx7yVFPO8BvwCbjM9z7DPXI6ic0wtYZny8DPif+UJ5SEq5G9XKJKW0Yu0FLJfKAaCYEKJcrgSaijRiT0svYJWUMlZKGQwEAk1MFlw6pJShUsqjxsf3gLNAeSz8c08n7rRY0mcupZRRxqe2xk0CHYA1xtcf/8yT/lusAZ4TIqmtde5KJ/a0WMS/FwAhRAWgG/Cj8bkgBz9znaCyRwJ/CSGOCNWOHqCslDLU+Pg6UNY8oWVKWrGWBy6n2O8K6f+CMpdRxqmNn1JMpVpk7MZpjAaov4rzzOf+WNyQBz5z41TTceAmsAM1oouQUiYYd0kZX3LsxvfvAiVzNeAUHo9dSpn0uX9u/NxnCyHsja9Z0uc+BxgHGIzPS5KDn7lOUNnTSkrZEOgKvC2EaJPyTanGsHli/X5eitXoe6A64AGEAl+ZNZp0CCGcgT+AsVLKyJTvWfLnnkrceeIzl1ImSik9gAqokVxt80aUeY/HLoSoB0xA/QyNgRLAePNF+CQhRHfgppTyiKnOoRNUNkgprxq/3gTWof7PcCNpmG38etN8EWYorVivAhVT7FfB+JrFkFLeMP6f2QAs4uGUkkXFLoSwRf2SXymlXGt82eI/99TiziufeRIpZQTgAzRHTX8l1RxNGV9y7Mb3iwLhuRvpk1LE3sU45SqllLHAEizvc28J9BRChACrUFN735CDn7lOUFkkhHASQhROegx0Bk4BG4Ghxt2GAhvME2GmpBXrRuAV4yqhZsDdFFNSFuGxufbeqM8eVOwDjCuFqgI1gUO5HR8kz8MvBs5KKb9O8ZZFf+5pxZ1HPvPSQohixseOQCfUNTQfoK9xt8c/86T/Fn2BncZRba5LI/b/UvwxI1DXcVJ+7mb/9yKlnCClrCClrIJa9LBTSjmInPzMc2ulR37ZgGqolUsngNPAJ8bXSwL/AOeBv4ES5o7VGNevqGmZeNR88OtpxYpaFTQfNXd/EvC0wNhXGGPzN/6DL5di/0+MsZ8Dupox7lao6Tt/4Lhxe8HSP/d04s4Ln7kbcMwY4yngU+Pr1VBJMxD4HbA3vu5gfB5ofL+aBca+0/i5nwJ+5uFKP4v49/LYz9COh6v4cuwz16WONE3TNIukp/g0TdM0i6QTlKZpmmaRdILSNE3TLJJOUJqmaZpF0glK0zRNs0g6QWkFihCiZIrq0NfFo1W67cwdX0pCiHZCiBa5dK4QIUQp4+NGxiraDXLj3JqWFpuMd9G0/ENKGY4q2YMQwhuIklLOMlc8Qggb+bBu2ePaAVHAvhw6Xma+3w1VyLO/lPJYdo+jaTlBj6C0As84YvA1Fv/dnuIO/l3GIp1+QoizQojGQoi1QvVzmmrcp4oQ4j8hxErjPmuEEIUycdw5QvUSGyOE6CFUf5xjQoi/hRBljcVahwPvGkd3rYXqD9Q3RdxRxq/thBB7hBAbgTPGwqMzhRCHjYVG38rkR1EHWA8MkVKapSKEpqWkE5RW0AlgLqrvTiPgJ+DzFO/HSSk9UX1tNgBvA/UALyFEUiXmZ4HvpJR1gEhgpLGmXXrHtZNSekopvwL2As2klA1QNc3GSSlDjOecLaX0kFLuyeDnaAiMkVLWQlXcuCulbIwqNDrMWIooIxuAUVLKvZnYV9NMTk/xaQWdPSrh7FAlz7BGlVdKstH49SRwWhprngkhLqAKX0YAl6WU/xr3+xl4B9iWwXFXp3hcAVhtHGHZAcHZ+DkOSdWTCVR9SLcUo62iqDp5GR33b+ANIcR2KWViNmLQtBylE5RW0AlU4mmexvuxxq+GFI+Tnif9/+fxemEyE8e9n+LxXOBrKeVGIUQ7VAfb1CRgnPUQQlihkllqxxPAaCnl9jSOk5ZRqFHbd0BmpwU1zWT0FJ9W0MUCpYUQzUG1mxBCuGbxGJWSvh94GTVldy4Lxy3Kw5YEQ1O8fg/Vej1JCNDI+LgnqvNqarYDI4zTjAghahkr7yOE+C+dn8NgjL+2EGJKOvtpWq7QCUor6Ayo0v9fCiFOoCp4Z3Vp9zlU48qzQHHgeyllXBaO6w38LoQ4AtxK8fqfQO+kRRKoXkxtjcdrzqOjppR+BM4AR4UQp4AfABvjMvJ0W2xLKWNQya+nEOLt9H9sTTMtXc1c056CcbXdJillPXPHkhGhOqBWk1J+a+5YNC0z9DUoTSsgpJSbzB2DpmWFHkFpmqZpFklfg9I0TdMskk5QmqZpmkXSCUrTNE2zSDpBaZqmaRZJJyhN0zTNIv0/BXDIjyB5EIYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "HeatCapacity = constant_volume_heat_capacity\n",
    "\n",
    "R = const.N_A * const.KB_J_PER_K\n",
    "plt.plot(temps, HeatCapacity(eq, temps) / R, \"-k\", label=\"Cv_rot Equilibrium\")\n",
    "plt.plot(temps, HeatCapacity(para, temps) / R, \"-b\", label=\"Cv_rot Para\")\n",
    "plt.plot(temps, HeatCapacity(ortho, temps) / R, \"-r\", label=\"Cv_rot Ortho\")\n",
    "plt.plot(\n",
    "    temps,\n",
    "    0.25 * HeatCapacity(para, temps) / R + 0.75 * HeatCapacity(ortho, temps) / R,\n",
    "    \"-g\",\n",
    "    label=\"Cv_rot 1:3 para:ortho\",\n",
    ")\n",
    "plt.plot(df_brink_T, df_brink_Cv, \"+g\")\n",
    "plt.plot(df_eucken_T, df_eucken_Cv, \"+g\")\n",
    "plt.plot(df_gia_T, df_gia_Cv, \"+g\")\n",
    "plt.plot(df_parting_T, df_parting_Cv, \"+g\")\n",
    "plt.plot(df_ce_T, df_ce_Cv, \"+g\", label=\"experimental data\")\n",
    "plt.legend(loc=\"upper right\", frameon=False)\n",
    "plt.xlim(10, 400)\n",
    "plt.ylim(-0.1, 2.8)\n",
    "plt.xlabel(\"Temperature, K\")\n",
    "plt.ylabel(\"Cv (rotational)/R\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>None</td></tr><tr><td>Terra</td><td>0.17.0.dev0+7dde867</td></tr><tr><td>Aer</td><td>0.8.0</td></tr><tr><td>Ignis</td><td>0.6.0.dev0+1537c75</td></tr><tr><td>Aqua</td><td>0.9.0.dev0+f7e1440</td></tr><tr><td>IBM Q Provider</td><td>0.12.0.dev0+8f3168b</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.7.9 (default, Aug 31 2020, 07:22:35) \n",
       "[Clang 10.0.0 ]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>2</td></tr><tr><td>Memory (Gb)</td><td>12.0</td></tr><tr><td colspan='2'>Thu Apr 01 09:50:07 2021 EDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2021.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
